text:
% source code to produce analytical predictions
% for 2d fast variance predictions paper for IEEE.
%
% Empirical and FFT_based predictions are computed and saved
%
% Copyright by Yingying Zhang-O'Connor, University of Michigan 2006

clear
clc

% -------------------------------
% set up parameters
% -------------------------------
printm 'set up parameters'
% calibration for single integral
% true to compute on the fly; false to load the saved table
t.calibrate = true;

t.ellipse = false;
t.zubal = false;
t.ncat = true;
t.fov = 'ellipse'; % 'ellipse' or 'disc'

t.const = false; % standard: \til{R}_phi = kappa_{center}
t.kappa = true; % certainty-based \til{R}_phi = kappa_j

if t.zubal
	img = 1;
	nx = 128/img;
	ny = nx;
	nb = ceil(888/4/img);
	na = ceil(984/4/img);
	f.ray_spacing = 0.9765625*4*img; %mm
	f.strip_width = 0.9765625*4*img;
	f.pixel_size = 0.9765625*4*img;
elseif t.ncat
	img = 2;
	nx = 512/img;
	ny = nx;
	nb = ceil(888/img);
	na = ceil(984/img);
	f.ray_spacing = 0.9765625*img; %mm
	f.strip_width = 0.9765625*img;
	f.pixel_size = 0.9765625*img;
	f.channel_offset = 0;
end

% fan-beam parameters
f.orbit = 360;
f.orbit_start = 0;
f.dis_src_det = 949.075;
f.dis_iso_det = 408.075;

% phantom parameters
if t.ellipse
	m.water = 0.2; % cm^2/g
	m.bone = 0.1;
	d.ellipse = 0.7*m.water/10; % g's g/cc; pixel size is mm, so convert into mm^(-1)
	d.bone = (2.5-1)*m.bone/10;
	d.soft = 1*m.water/10;
	c_x = (nx+1)/2;
	c_y = (ny+1)/2;
	bx_center = 0;
	by_center = 0;
	bx_radius = (nx - 38)/2;
	by_radius = 25;
	disc_radius = 10;
end

% ------------------------------
% generate the true phantom
% ------------------------------
printm 'set up xtrue'
xtrue = zeros(nx,ny);
if t.ellipse
	x_center = [bx_center; bx_center+20];
	y_center = [by_center; by_center];
	x_radius = [bx_radius; disc_radius];
	y_radius = [by_radius; disc_radius];
	ampl = [d.ellipse; d.bone];
	params = [x_center y_center x_radius y_radius zeros(2,1) ampl];
	params(:,1:4) = params(:,1:4)*f.pixel_size;
	xtrue = zeros(nx,ny);
	xtrue = ellipses(nx, ny, params, f.pixel_size, f.pixel_size);

elseif t.zubal
	dir = '~fessler/l/src/matlab/alg/data';
	xtrue = read_zubal_attn(dir) * 2.1; % mulist = [0.002 0.0096 0.0120] *2.1 to make CT value;

	if nx < 128
		xtrue = downsample2(xtrue,img);
	end

elseif t.ncat
%	xtrue = fld_read('~yzz/data/phantoms/ncat/ncat,256,slice,140.fld')/2;
%	ddir = path_find_dir([filesep 'transmission']);
%	ddir = strrep(ddir, 'transmission', 'data');
%	xtrue = fld_read([ddir '/ncat,256,slice,140,ct,x100.fld']) / 2;
	ddir = '~fessler/l/dat/phantom,digital/unc/ncat1.1';
	xtrue = fld_read([ddir '/ncat,256,slice,140.fld']) / 2;

	if (nx < 1024) && size(xtrue,2) == 1024
		img = 1024/nx;
		xtrue = downsample2(xtrue,img);
	end
end

% -------------------------------
% generate mask
% -------------------------------
printm 'generate support'
if t.ellipse
	x_mask = bx_radius+5;
	y_mask = by_radius+5;
elseif t.zubal || t.ncat
	if streq(t.fov, 'ellipse')
		printm('ellipse FOV')
		x_mask = (nx - nx/16 + 10)/2;
		y_mask = (ny - (ny/4 + nx/16) + 50)/2;
	elseif streq(t.fov, 'disc')
		printm('circular FOV')
		ss = ((-(nb-1)/2:(nb-1)/2).' - f.channel_offset) * f.ray_spacing;
		smax = max(abs(ss));
		rmax = sin(smax/f.dis_src_det)*(f.dis_src_det - f.dis_iso_det);
		x_mask = floor((rmax-2)/2);
		y_mask = x_mask;
	end
end
f.support = sprintf('ellipse 0 0 %g %g', x_mask, y_mask);

% new code by jf showing how to generate system matrix
sg = sino_geom('fan', ...
	'nb', nb, 'na', na, ...
	'dsd', f.dis_src_det, ...
	'dod', f.dis_iso_det, ...
	'ds', f.ray_spacing, ...
	'strip_width', f.strip_width, ...
	'orbit', f.orbit, ...
	'orbit_start', f.orbit_start, ...
	'offset_s', 0);

ig = image_geom('nx', nx, 'ny', ny, 'dx', f.pixel_size);
ig.mask = ig.circ(max(x_mask, y_mask));
if x_mask ~= y_mask
	warn 'need to put in ellipse using ellipse_im'
end

G = Gtomo2_wtmex(sg, ig);

if 0 % this is an old way - do not use it
	sys = arg_pair('system', 14, ...
	'nx', nx, 'ny', ny, ...
	'nb', nb, 'na', na, ...
	'support',	f.support, ...
	'pixel_size',	f.pixel_size, ...
	'ray_spacing',	f.ray_spacing, ...
	'strip_width',	f.strip_width, ...
	'orbit',	f.orbit, ...
	'orbit_start',	f.orbit_start, ...
	'src_det_dis',	f.dis_src_det, ...
	'obj2det_x',	f.dis_iso_det, ...
	'obj2det_y',	f.dis_iso_det, ...
	'flip_y', 1, ...
	'scale', 0, ...
	'channel_offset', f.channel_offset, ...
	'source_offset', 0);
	Gdsc = Gtomo2_dscmex(sys, 'nthread', 2); % use it to generate the mask
	mask = Gdsc.arg.mask;
end

% -------------------------------
% system matrix
% -------------------------------
printm 'set up fan-beam system using DD projector'

if 0 % this is an old way - do not use it
	G = Gtomo2_dd(mask, nb, na, ...
	'pixel_size', f.pixel_size, ...
	'ray_spacing', f.ray_spacing, ...
	'strip_width', f.strip_width, ...
	'dis_src_det', f.dis_src_det, ...
	'dis_iso_det', f.dis_iso_det, ...
	'orbit', f.orbit, ...
	'orbit_start', f.orbit_start, ...
	'xscale', 1, ...
	'yscale', 1, ...
	'channel_offset', 0);
end

% ---------------------------------
% compute statistical weighting wi
% ---------------------------------
printm 'compute statistical weighting wi'
lb = reshape(G*xtrue(mask(:)), [nb na]);
bi = 10^6;
yb = bi*exp(-lb);
yi = poissrnd(yb);
wi = yb; % this is the true case. In reality no yb. Use yi plug-in method
lhat = -log(yi./bi);

W = diag_sp(wi(:));

% ---------------------------------
% choose \alpha for all penalty
% ---------------------------------
% regularization parameter
if nx == 128
	f.l2b = 12; % for 128x128
elseif nx == 32
	f.l2b = 14;
elseif (nx == 512) || (nx == 256)
	f.l2b = 11;
else
	f.l2b = 17;
end
f.beta = 2^f.l2b;

% ----------------------------------------
% Compute the H for all penalty cases
% ----------------------------------------
dx = f.pixel_size;
dy = dx;
np = G.arg.np; % number of pixels inside mask
ds = f.ray_spacing;
dp = 1; % distance power
dbeta = 2*pi/na;

printm('Polar frequency coordinates')
% decide rhomax
rhomax = 1 / (2 * dx);
drho = 1/nx/dx;
rho1 = (-nx/2:nx/2-1) / (nx * dx);
rho = rho1(nx/2 + 1:end);
rho = rho(:);
nrho = length(rho);
nPhi = na/4;
Phi = ((0:nPhi-1)')/nPhi * 2 * pi;
dPhi = abs(Phi(2)-Phi(1));

printm('compute w_0(Phi): angular dependent weighting')
nw = na/4;
phiw = ((0:nw-1)')/nw * pi; % to avoid = pi/2 since tan(pi/2) = inf
dphiw = abs(phiw(2)-phiw(1));
% w_0(Phi)
[wknot, tw] = adw_fan(G.arg, wi, mask, Phi);
sprintf('w_0(Phi) pre-compute: %g seconds', tw)

%
% pre-compute G(\rho, \Phi) in H(\rho, \Phi) at image center
%
% Jacobian: J(sj) \approx J(s_center) at image center: Ds0/Dsd
Jj = G.arg.dis_src_iso/G.arg.dis_src_det;

if streq(t.fov, 'ellipse')
	if 1 % d_0(\phi) at image center
		fov = x_mask .* y_mask ./ sqrt((x_mask .* cos(phiw)).^2 ...
			+ (y_mask .* sin(phiw)).^2);
		fov = 2 .* fov(:) .* dx;
	end

elseif streq(t.fov, 'disc')
	fov = sum(mask(:,nx/2)) .* dx;
	fov = repmat(fov, [length(phiw) 1]);
end

tt = cputime;
G_center = zeros(nrho, nPhi);
for ia = 1:length(phiw)
	G_center = G_center + fov(ia).* sinc(rho * (fov(ia) .* sin(Phi-phiw(ia))).').^2;
end
sprintf('G_0(rho, Phi) at center pre-compute: %g seconds', cputime - tt)
% A_0(\rho,\Phi) at image center
A = dx .* sinc(dx .* (rho * cos(Phi).')) .* ...
	dy .* sinc(dx .* (rho * sin(Phi).')) .* ...
	repmat((sinc(ds .* Jj .* rho)), [1 length(Phi)]); % Nrho x Phi

Ga = G_center .* (A.^2) .* dphiw ./ abs(ds .* dx .* dy .* dbeta);

% -----------------------------------------
% Compute the approx standard deviation map
% -----------------------------------------
if t.const || t.kappa

	kappan = G' * wi(:);
	kappad = G' * (ones(size(wi(:))));
	kappa = sqrt( kappan ./ kappad );

	f.data_dir = 'data/standard_penal/256image/';

	if t.const
	printm('standard quadratic penalty')
	printm('kappa is a const for all pixels')
	kappa_r = embed(kappa, mask);
	kappa_cen = kappa_r(nx/2+1, ny/2+1) .* ones(size(kappa));
	kappa_v = kappa; % varies for pixel

	his = f.l2b + log2(kappa_v.^2) - log2(kappa_cen.^2);
	his = hist(his, (6:0.5:14).'); % histogram
	clear kappa_r

	kappa = kappa_cen;
	kappa_r = embed(kappa, mask);

	% empirical and fft_based for comparison
	f.iot = [f.data_dir 'std_iot,iter71,ntime250,l2b11,256,nblock41,const.fld'];
	std_iot = fld_read([f.data_dir 'std_iot,iter71,ntime250,l2b11,256,nblock41,const.fld']);
	std_fft = fld_read([f.data_dir 'std_vh,dbsized_mask,iter71,ntime250,l2b11,256,nblock41,const.fld']);

	elseif t.kappa
	printm('certainty-base quadratic penalty')
	printm('kappa varies for pixels')
	kappa_r = embed(kappa, mask);

	% empirical and fft_based for comparison
%	std_iot = fld_read([f.data_dir 'std_iot,iter71,ntime250,l2b11,256,nblock41,kappa.fld']);
	std_iot = zeros(size(kappa_r)); % fake jf
%	std_fft = fld_read([f.data_dir 'std_vh,dbsized_mask,iter71,ntime250,l2b11,256,nblock41,kappa.fld']);
	std_fft = zeros(size(kappa_r)); % fake jf
	end

	r_sum_con = 1 + sqrt(2).^(2-dp);
	nl = [1 0 1 1]; ml = [0 -1 -1 1];
	phil = atan2(ml,nl);

	uu = rho * cos(Phi).';
	vv = rho * sin(Phi).';
	R1 = 0;
	for il = 1:4
	Rl = 4 .* sin((uu .* nl(il) + vv .* ml(il))*2*pi*dx/2).^2 ./ (sqrt(nl(il)^2 + ml(il)^2)).^dp;
	R1 = R1 + Rl; % rho x Phi
	end

	printm('double summation')
	tt = cputime;
	for ia = 1:size(wknot,2)
	ticker(mfilename, ia, size(wknot,2))
	Hd = wknot(:,ia) * Ga(:,ia).'; % Np x Nrho
	Hu = wknot(:,ia) * (Ga(:,ia) .* rho).'; % Np x Nrho

	R = kappa.^2 * R1(:,ia).'; % Np x Nrho
	int_rho(:,ia) = rhomax .* mean(Hu ./ ((Hd + f.beta .* R).^2), 2);
	end
	var_pre = mean(int_rho,2) .* (2*pi) .* (dx .* dy);
	std_pre = embed(sqrt(var_pre), mask);
	sprintf('db-integral prediction: %g second', cputime - tt)

	if t.calibrate
	% use extended mask for FFT in calibration
	maskcal = true([nx ny]);
	Gcal = Gtomo2_dd(maskcal, nb, na, ...
	'pixel_size', f.pixel_size, ...
	'ray_spacing', f.ray_spacing, ...
	'strip_width', f.strip_width, ...
	'dis_src_det', f.dis_src_det, ...
	'dis_iso_det', f.dis_iso_det, ...
	'orbit', f.orbit, ...
	'orbit_start', f.orbit_start, ...
	'xscale', 1, ...
	'yscale', 1, ...
	'channel_offset', 0);

	c.l2b = [2:0.15:20].'; % column
	c.beta = 2.^c.l2b;
	for indx = 1:length(c.beta)
	c.R0 = Robject(maskcal, 'edge_type', 'tight', 'beta', c.beta(indx), ...
	'type_denom', 'matlab', 'potential', 'quad', 'distance_power', dp);
	% mask size matters for smaller \beta
	[psf0, c.var0(indx), fwhm_t(indx)] = qpwls_psf(Gcal, c.R0, 1, maskcal, [], [0 0]);
	end
	c.var0 = c.var0(:); % column
	fwhm_t = fwhm_t(:);

	ej = zeros([nx ny]);
	ej(nx/2+1,ny/2+1) = 1;
	center_in = find(ej(mask(:))==1);
	[wknotcal, tw] = adw_fan(G.arg, ones(size(wi)), mask, Phi);
	CC = (dx .* dy) ./ (ds);
	KK = CC/(rhomax)^3;
	for indx = 1:length(c.beta)
	c.var_pre(indx) = mean(((dx .* dy) .* (2*pi) ./ 3) ./ ...
	(KK .* wknotcal(center_in,:) .* 2 + (2*pi*dx)^2*c.beta(indx)*r_sum_con));
	end
	c.var_pre = c.var_pre(:);

	ratio = sqrt(c.var0(:) ./ c.var_pre);

%	fld_write('data/cal,var_fft,fwhm,qpuls,l2b2to20,inl2b,2w0.fld', [c.var0 fwhm_t]);
%	fld_write('data/cal,ratio,in_std,qpuls,l2b2to20,inl2b,2w0.fld', ratio)
	else
		% load the calibration scale factor needed for single integral
		% Calibration is done for a range of fwhms
		% about (1.13)^2 in our case
		tmp = fld_read('data/cal,var_fft,fwhm,qpuls,l2b2to20,inl2b,2w0.fld');
		fwhm_t = tmp(:,2);
		ratio = fld_read('data/cal,ratio,in_std,qpuls,l2b2to20,inl2b,2w0.fld');
	end
	Rr = Robject(kappa_r, 'edge_type', 'tight', 'beta', f.beta, ...
		'type_denom', 'matlab', 'potential', 'quad', 'distance_power', dp);
	[psf, var, fwhm] = qpwls_psf(G, Rr, 1, mask, W, [0 0]);
	f_loc = find(abs(fwhm - fwhm_t) == min(abs(fwhm - fwhm_t)));

	CC = (dx .* dy) ./ (ds);
	KK = CC/(rhomax)^3;
	tt = cputime;
	var_pre1 = mean(1 ./ (KK .* wknot * 2 + (2*pi*dx)^2*f.beta*r_sum_con .* repmat(kappa.^2,[1 length(Phi)])), 2);
	std_pre1 = embed(sqrt(var_pre1 .* (dx .* dy) .* (2*pi) ./ 3), mask);
	std_pre1 = std_pre1 .* ratio(f_loc);
	sprintf('one-integral prediction: %g second', cputime - tt)
end

% jf: display results
im pl 2 2
im(1, std_pre), cbar
im(2, std_pre1), cbar

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%--- help for load_data ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named load_data
%
%       dsge/load_data    generic/load_data
%

author: jmaih


text:
function r = atanh(a)
%ATANH        Hessian (elementwise) inverse hyperbolic tangent
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = atanh(a.x);
    f = 1 / ( 1 - sqr(a.x) );
    r.dx = a.dx * f;
    r.hx = ( a.hx + reshape( (a.x*r.dx) * a.dx.' , size(a.hx) ) ) * f;

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = atanh(a.x);
    if issparse(a.hx)               % input sparse

      ax = 1 ./ ( 1 - sqr(full(a.x(:))) );
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = a.x(:).*ax;
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,r.dx);
      end
      [ia,ja,sa] = find(a.hx);      % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = atanh(a.x);
      ax = a.x(:).';
      f = 1 ./ ( 1 - sqr(ax) );
      f = f(ones(N*N,1),:);
      r.dx = a.dx .* f(1:N,:);
      adx = repmat(ax,N,1) .* r.dx;
      r.hx = ( a.hx + adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:) ) .* f;

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [rfrac,isonfly] = center_of_rotation2(trx,fly,debug)
if ~exist('debug','var')
  debug = false;
end
N = 100;

% note that if dtheta = 0, this will produce 0

cost = cos(trx(fly).theta_mm);
sint = sin(trx(fly).theta_mm);
dacost = 2*diff(trx(fly).a_mm.*cost);
dbcost = 2*diff(trx(fly).b_mm.*cost);
dasint = 2*diff(trx(fly).a_mm.*sint);
dbsint = 2*diff(trx(fly).b_mm.*sint);
Z = dacost .* dbcost + dbsint.*dasint;
Minv = zeros(2,2,trx(fly).nframes-1);
Minv(1,1,:) = dbcost ./ Z;
Minv(1,2,:) = dbsint ./ Z;
Minv(2,1,:) = -dasint ./ Z;
Minv(2,2,:) = dacost ./ Z;
Minv = reshape(Minv,[4,trx(fly).nframes-1]);
x = [ diff(trx(fly).x_mm) ; diff(trx(fly).y_mm) ];
rfrac = -[Minv(1,:) .* x(1,:) + Minv(3,:) .* x(2,:);...
  Minv(2,:) .* x(1,:) + Minv(4,:) .* x(2,:)];
% when no rotation, set center of rotation to middle of fly by default
rfrac(isnan(rfrac)) = 0;
rfrac0 = rfrac;
isoutofbounds = sum(rfrac.^2,1) > 1;
psi = linspace(0,2*pi,N)';
cospsi = cos(psi);
sinpsi = sin(psi);
idx = find(isoutofbounds);
nout = length(idx);
if nout > 0,
  x1 = repmat(trx(fly).x_mm(idx),[N,1]) + repmat(2*trx(fly).a_mm(idx).*cost(idx),[N,1]).*repmat(cospsi,[1,nout]) - ...
    repmat(2*trx(fly).b_mm(idx).*sint(idx),[N,1]).*repmat(sinpsi,[1,nout]);
  y1 = repmat(trx(fly).y_mm(idx),[N,1]) + repmat(2*trx(fly).a_mm(idx).*sint(idx),[N,1]).*repmat(cospsi,[1,nout]) + ...
    repmat(2*trx(fly).b_mm(idx).*cost(idx),[N,1]).*repmat(sinpsi,[1,nout]);
  x2 = repmat(trx(fly).x_mm(idx+1),[N,1]) + repmat(2*trx(fly).a_mm(idx+1).*cost(idx+1),[N,1]).*repmat(cospsi,[1,nout]) - ...
    repmat(2*trx(fly).b_mm(idx+1).*sint(idx+1),[N,1]).*repmat(sinpsi,[1,nout]);
  y2 = repmat(trx(fly).y_mm(idx+1),[N,1]) + repmat(2*trx(fly).a_mm(idx+1).*sint(idx+1),[N,1]).*repmat(cospsi,[1,nout]) + ...
    repmat(2*trx(fly).b_mm(idx+1).*cost(idx+1),[N,1]).*repmat(sinpsi,[1,nout]);
  d = (x1 - x2).^2 + (y1 - y2).^2;
  [~,j] = min(d,[],1);
  rfrac(1,idx) = cospsi(j);
  rfrac(2,idx) = sinpsi(j);
end
isonfly = ~isoutofbounds;

if debug,
  ntry = 50;
  psitry = linspace(0,2*pi,ntry+2);
  psitry = psitry(2:end-1);
  [rhotry,psitry] = meshgrid(linspace(0,1,ntry),psitry);
  rtry = [rhotry(:).*cos(psitry(:)),rhotry(:).*sin(psitry(:))]';
  rtry = [rtry,[0;0]];
  colorsplot = 'rg';
  colorsplot2 = 'mc';
  [~,order] = sort(-abs(modrange(diff(trx(fly).theta_mm),-pi,pi)));
  for i = order,
    u = [0,0]; v = [0,0];
    uout = [0,0]; vout = [0,0];
    utry = zeros(2,length(rtry)); vtry = zeros(2,length(rtry));
    clf; hold on;
    for j = [i,i+1],
      u(j-i+1) = trx(fly).x_mm(j) + rfrac(1,i)*trx(fly).a_mm(j)*2*cos(trx(fly).theta_mm(j)) - rfrac(2,i)*trx(fly).b_mm(j)*2*sin(trx(fly).theta_mm(j));
      v(j-i+1) = trx(fly).y_mm(j) + rfrac(1,i)*trx(fly).a_mm(j)*2*sin(trx(fly).theta_mm(j)) + rfrac(2,i)*trx(fly).b_mm(j)*2*cos(trx(fly).theta_mm(j));
      uout(j-i+1) = trx(fly).x_mm(j) + rfrac0(1,i)*trx(fly).a_mm(j)*2*cos(trx(fly).theta_mm(j)) - rfrac0(2,i)*trx(fly).b_mm(j)*2*sin(trx(fly).theta_mm(j));
      vout(j-i+1) = trx(fly).y_mm(j) + rfrac0(1,i)*trx(fly).a_mm(j)*2*sin(trx(fly).theta_mm(j)) + rfrac0(2,i)*trx(fly).b_mm(j)*2*cos(trx(fly).theta_mm(j));
      utry(j-i+1,:) = trx(fly).x_mm(j) + rtry(1,:)*trx(fly).a_mm(j)*2*cos(trx(fly).theta_mm(j)) - rtry(2,:)*trx(fly).b_mm(j)*2*sin(trx(fly).theta_mm(j));
      vtry(j-i+1,:) = trx(fly).y_mm(j) + rtry(1,:)*trx(fly).a_mm(j)*2*sin(trx(fly).theta_mm(j)) + rtry(2,:)*trx(fly).b_mm(j)*2*cos(trx(fly).theta_mm(j));
      ellipsedraw(trx(fly).a_mm(j)*2,trx(fly).b_mm(j)*2,trx(fly).x_mm(j),trx(fly).y_mm(j),trx(fly).theta_mm(j),colorsplot(j-i+1));
      plot(trx(fly).x_mm(j)+trx(fly).a_mm(j)*2*cos(trx(fly).theta_mm(j)),trx(fly).y_mm(j)+trx(fly).a_mm(j)*2*sin(trx(fly).theta_mm(j)),'x','color',colorsplot(j-i+1));
      htest(j-i+1) = plot([trx(fly).x_mm(j),u(j-i+1),uout(j-i+1)],[trx(fly).y_mm(j),v(j-i+1),vout(j-i+1)],'o-','color',colorsplot(j-i+1),'markerfacecolor',colorsplot(j-i+1)); %#ok<AGROW>
    end
    d = (utry(1,:)-utry(2,:)).^2 + (vtry(1,:)-vtry(2,:)).^2;
    [dexp,k] = min(d);
    uexp = utry(:,k);
    vexp = vtry(:,k);
    dexp = sqrt(dexp);
    dtest = sqrt(diff(u).^2 + diff(v).^2);
    dout = sqrt(diff(uout).^2 + diff(vout).^2);
    plot(uexp,vexp,'k.-');
    plot(u,v,'b.-');
    plot(uout,vout,'b.-');
    for j = [i,i+1],
      hexp(j-i+1) = plot([uexp(j-i+1),trx(fly).x_mm(j)],[vexp(j-i+1),trx(fly).y_mm(j)],'o-','color',colorsplot2(j-i+1),'markerfacecolor',colorsplot2(j-i+1)); %#ok<AGROW>
    end
    ax = nan(1,4);
    ax1 = nan(1,4);
    for j = [i,i+1],
      [ax1(1),ax1(2),ax1(3),ax1(4)] = ellipse_to_bounding_box(trx(fly).x_mm(j),trx(fly).y_mm(j),2*trx(fly).a_mm(j),2*trx(fly).b_mm(j),trx(fly).theta_mm(j));
      ax([1,3]) = min(ax([1,3]),ax1([1,3]));
      ax([2,4]) = max(ax([2,4]),ax1([2,4]));
    end
    ax = ax + [-3,3,-3,3];
    legend([htest,hexp],'analytic1','analytic2','empirical1','empirical2');
    title(sprintf('red = frame %d, green = frame %d, dtest = %f, dout = %f, dexp = %f',i,i+1,dtest,dout,dexp));
    axis equal
    axis(ax)
    input('');
  end
end

author: kristinbranson

repo: JAABA


text:
function delay = fast_slicedelay(TR,nSlices,Slice,AcqOrder)
% delay = fast_slicedelay(TR,nSlices,Slice,AcqOrder)
%
% Computes the amount of time after the start of a TR that a slice
% is aquired given the TR, the number of slices, and the order of
% acquisition (linear or interleaved).
%

if(nargin ~= 4)
     msg = 'USAGE: delay = fast_slicedelay(TR,nSlices,Slice,AcqOrder)';
  qoe(msg);error(msg);
end

if( ~strcmpi(AcqOrder,'linear') & ~strcmpi(AcqOrder,'interleaved'))
  msg = sprintf('AcqOrder = %s, must be either linear or interleaved',AcqOrder);
  qoe(msg);error(msg);
end

dt = TR/nSlices;

if( strcmpi(AcqOrder,'linear') )
  SliceOrder = [0:nSlices-1];
end

if( strcmpi(AcqOrder,'interleaved') )
  SliceOrder = [[0:2:nSlices-1] [1:2:nSlices] ];
end

nthAcq = find(SliceOrder == Slice) - 1;
delay = dt*nthAcq;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [msg, niter] = parallel_protocol(engine, evidence, msg)

bnet = bnet_from_engine(engine);
N = length(bnet.dag);
ns = bnet.node_sizes(:);

if ~isempty(engine.filename)
  fid = fopen(engine.filename, 'w');
  if fid == 0
    error(['could not open ' engine.filename ' for writing'])
  end
else
  fid = [];
end

converged = 0;
iter = 1;
hidden = find(isemptycell(evidence));
bel = cell(1,N);
old_bel = cell(1,N);
%nodes = mysetdiff(1:N, engine.disconnected_nodes);
nodes = find(~engine.disconnected_nodes_bitv);
while ~converged & (iter <= engine.max_iter)
  % Everybody updates their state in parallel
  for n=nodes(:)'
    cs_msg = children(engine.msg_dag, n);
    %msg{n}.lambda = compute_lambda(n, cs, msg);
    msg{n}.lambda = prod_lambda_msgs(n, cs_msg, msg, engine.msg_type);
    ps_orig = parents(bnet.dag, n);
    msg{n}.pi = CPD_to_pi(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, msg, evidence);
  end

  changed = 0;
  if ~isempty(fid)
    fprintf(fid, 'ITERATION %d\n', iter);
  end
  for n=hidden(:)' % this will not contain any disconnected nodes
    old_bel{n} = bel{n};
    bel{n}  = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda);
    if ~isempty(fid)
      fprintf(fid, 'node %d: %s\n', n, bel_to_str(bel{n}, engine.msg_type));
    end
    if engine.storebel
      engine.bel{n,iter} = bel{n};
    end
    if (iter == 1) | ~approxeq_bel(bel{n}, old_bel{n}, engine.tol, engine.msg_type)
      changed = 1;
    end
  end
  %converged = ~changed;
  converged = ~changed & (iter > 1);  % Sonia Leach changed this

  if ~converged
    % Everybody sends to all their neighbors in parallel
    for n=nodes(:)'
      % lambda msgs to parents
      ps_msg = parents(engine.msg_dag, n);
      ps_orig = parents(bnet.dag, n);
      for p=ps_msg(:)'
	j = engine.child_index{p}(n); % n is p's j'th child
	old_msg = msg{p}.lambda_from_child{j}(:);
	new_msg = CPD_to_lambda_msg(bnet.CPD{bnet.equiv_class(n)}, engine.msg_type, n, ps_orig, ...
				    msg, p, evidence);
	lam_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type);
	msg{p}.lambda_from_child{j} = lam_msg;
      end

      % pi msgs to children
      cs_msg = children(engine.msg_dag, n);
      for c=cs_msg(:)'
	j = engine.parent_index{c}(n); % n is c's j'th parent
	old_msg = msg{c}.pi_from_parent{j}(:);
	%new_msg = compute_pi_msg(n, cs, msg, c));
	new_msg = compute_bel(engine.msg_type, msg{n}.pi, prod_lambda_msgs(n, cs_msg, msg, engine.msg_type, c));
	pi_msg = convex_combination_msg(old_msg, new_msg, engine.momentum, engine.msg_type);
	msg{c}.pi_from_parent{j} = pi_msg;
      end
    end
    iter = iter + 1;
  end
end

if fid > 0, fclose(fid); end
%niter = iter - 1;
niter = iter;

%%%%%%%%%%

function str = bel_to_str(bel, type)

switch type
 case 'd', str = sprintf('%9.4f ', bel(:)');
 case 'g', str = sprintf('%9.4f ', bel.mu(:)');
end

%%%%%%%

function a = approxeq_bel(bel1, bel2, tol, type)

switch type
 case 'd', a = approxeq(bel1, bel2, tol);
 case 'g', a = approxeq(bel1.mu, bel2.mu, tol) & approxeq(bel1.Sigma, bel2.Sigma, tol);
end

%%%%%%%

function msg = convex_combination_msg(old_msg, new_msg, old_weight, type)

switch type
 case 'd', msg = old_weight * old_msg + (1-old_weight)*new_msg;
 case 'g', msg = new_msg;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
filename = 'CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
constraint_case = {'EQUALITY'};
optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';
incrementFactor = 1;
kfrac = 1.05;
designVariable = 'Density';
filterType = 'P1';

nsteps = 2;
Vfrac_final = 0.8;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
Perimeter_target = 5;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

% For all tests
plotting = false;
printing = false;
printing_physics = false;
monitoring = false;
maxiter = 5;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef PTKTreeSegment < PTKTree
    % PTKTreeSegment. A data structure representing a segmented airway tree
    %
    %     A root PTKTreeSegment is returned by
    %     PTKAirwayRegionGrowingWithExplosionControl. From this you
    %     can extract and analyse the resulting airway tree.
    %
    %     PTKTreeSegment is used in the construction and storage of
    %     the airway trees. A PTKTreeSegment stores an individual
    %     segment of the centreline tree, with references to the parent and child
    %     PTKTreeSegments, so that it is possible to reconstruct the entire
    %     tree from a single segment.
    %
    %     The way the airway voxels are stored in each segment is as follows:
    %         The wavefront only exists during the region growing. It is a thick
    %         layer of voxels which exists at the end of each segment
    %         which is currently growing. If the wavefront forms more than one
    %         connected component, the segment is ended and new child segments
    %         formed from the wavefront components. When new voxels are added to
    %         the front of the wavefront, old voxels from the back of the
    %         wavefront are pushed out into the Pending voxels of the segment.
    %         Pending voxels are 'accepted' or 'rejected' according to whether
    %         the heuristics have determined an explosion has occurred.
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        Colour         % A colourmap index allocated to this branch
    end

    properties (SetAccess = private)
        % Generation of this segment, starting at 1
        GenerationNumber

        % Indicates if this branch was terminated due to the generation number
        % being too high, which may indicate leakage
        ExceededMaximumNumberOfGenerations = false

        MarkedExplosion = false

    end

    properties (Access = private)

        % List of accepted (i.e. non-exploded) voxels
        AcceptedVoxelIndices

        % List of rejected (exploded) voxels
        RejectedVoxelIndices

        % Indices are allocated to the segment from the back of the wavefront
        % They are pending until the explosion control heuristic has determined
        % they are 'Accepted' or 'Exploded'.
        PendingVoxelIndices

        % Additional voxels that were not originally part of the segment but
        % were added later after a morphological closing operation
        ClosedPoints

        PreviousMinimumVoxels
        NumberOfVoxelsSkipped = 0
        LastNumberOfVoxels = 0

        PermittedVoxelSkips = 0
        IsFirstSegment

        % We never use less than this value at each step when computing the minimum wavefront
        % size over the image
        MinimumNumberOfPointsThreshold

        ExplosionMultiplier

    end

    methods
        function obj = PTKTreeSegment(parent, min_number_of_points_threshold, explosion_multiplier)
            if nargin > 0
                obj.Parent = parent;
                obj.MarkedExplosion = false;
                obj.PendingVoxelIndices = int32([]);
                obj.AcceptedVoxelIndices = int32([]);
                obj.RejectedVoxelIndices = int32([]);

                obj.ExplosionMultiplier = explosion_multiplier;
                obj.MinimumNumberOfPointsThreshold = min_number_of_points_threshold;

                if ~isempty(parent)
                    parent.AddChild(obj);
                    obj.PreviousMinimumVoxels = parent.PreviousMinimumVoxels;
                    obj.LastNumberOfVoxels = obj.PreviousMinimumVoxels;
                    obj.IsFirstSegment = false;
                    obj.GenerationNumber = parent.GenerationNumber + 1;
                else
                    obj.PreviousMinimumVoxels = 1000;
                    obj.LastNumberOfVoxels = obj.PreviousMinimumVoxels;
                    obj.IsFirstSegment = true;
                    obj.GenerationNumber = 1;
                end
            end
        end

        % Returns a list of voxels which have been accepted as not explosions.
        function accepted_voxels = GetAcceptedVoxels(obj)
            accepted_voxels = obj.ConcatenateVoxels(obj.AcceptedVoxelIndices);
        end

        % Returns the wavefront for this segment, which includes voxels that may
        % be separated into child segments
        function rejected_voxels = GetRejectedVoxels(obj)
            rejected_voxels = obj.ConcatenateVoxels(obj.RejectedVoxelIndices);
        end

        function endpoints = GetEndpoints(obj)
            endpoints = obj.AcceptedVoxelIndices{end};
        end

        % Returns all accepted region-growing points, plus those added from the airway closing operation
        function all_points = GetAllAirwayPoints(obj)
            all_points = [obj.GetAcceptedVoxels; obj.ClosedPoints];
        end

        % Points which are added later to close gaps in the airway tree
        function AddClosedPoints(obj, new_points)
            obj.ClosedPoints = [obj.ClosedPoints, new_points];
        end

        function RecomputeGenerations(obj, new_generation_number)
            obj.GenerationNumber = new_generation_number;
            children = obj.Children;
            for child = children
                child.RecomputeGenerations(new_generation_number + 1);
            end
        end

        % Returns the number of branches in this tree, from this branch
        % and excluding generations above the max_generation_number
        function number_of_branches = CountBranchesUpToGeneration(obj, max_generation_number)
            number_of_branches = 0;
            if obj.GenerationNumber > max_generation_number
                return;
            else
                branches_to_do = obj;
                while ~isempty(branches_to_do)
                    branch = branches_to_do(end);
                    branches_to_do(end) = [];
                    if branch.GenerationNumber <= max_generation_number
                        branches_to_do = [branches_to_do, branch.Children];
                        number_of_branches = number_of_branches + 1;
                    end
                end
            end
        end

        function AddColourValues(obj, new_colour_value)
            obj.Colour = new_colour_value;
            children = obj.Children;
            for child = children
                new_colour_value = new_colour_value + 1;
                if mod(new_colour_value, 7) == 0
                    new_colour_value = new_colour_value + 1;
                end
                child.AddColourValues(new_colour_value);
            end
        end

        function AddPendingVoxels(obj, indices_of_new_points)

            if isempty(obj.PendingVoxelIndices)
                obj.LastNumberOfVoxels = max(1, round(numel(indices_of_new_points)/2));
            end
            obj.PendingVoxelIndices{end + 1} = indices_of_new_points;

            if obj.MarkedExplosion
                obj.RejectAllPendingVoxelIndices;
                return;
            end

            number_of_points = numel(indices_of_new_points);

            if (number_of_points < obj.PreviousMinimumVoxels) && (~obj.IsFirstSegment)
                obj.PreviousMinimumVoxels = max(number_of_points, obj.MinimumNumberOfPointsThreshold);
            end

            if (number_of_points < obj.ExplosionMultiplier*obj.PreviousMinimumVoxels)
                obj.NumberOfVoxelsSkipped = 0;
            else
                obj.NumberOfVoxelsSkipped = obj.NumberOfVoxelsSkipped + 1;
            end

            % Keep track of the point at which an explosion starts to occur
            if (number_of_points <= obj.LastNumberOfVoxels)
                obj.LastNumberOfVoxels = number_of_points;
                obj.AcceptAllPendingVoxelIndices;
            end

            % Explosion control: we allow a certain number of consecutive
            % points to exceed the expansion limit.
            % Once exceeded, this segment is not permitted to expand further, and it is also marked so that it can be
            % deleted later.
            if (obj.NumberOfVoxelsSkipped > obj.PermittedVoxelSkips)
                obj.MarkedExplosion = true;
                obj.RejectAllPendingVoxelIndices;
            end
        end

        function EarlyTerminateBranch(obj)
            obj.CompleteThisSegment;
            if ~obj.MarkedExplosion
                obj.ExceededMaximumNumberOfGenerations = true;
            end
        end

        function CompleteThisSegment(obj)
            if obj.MarkedExplosion
                obj.RejectAllPendingVoxelIndices;
            else
                obj.AcceptAllPendingVoxelIndices;
            end
        end

        function MergeWithChild(obj)
            if (length(obj.Children) ~= 1)
                error('MergeWithChild should only be called if there is one child segment');
            end
            child = obj.Children(1);
            if ~isempty(child.PendingVoxelIndices)
                error('Programming error: MergeWithParent should not be called while a segment is in progress');
            end
            grandchildren = child.Children;
            obj.Children = [setdiff(obj.Children, child), grandchildren];
            for grandchild = obj.Children
                grandchild.Parent = obj;
                grandchild.RecomputeGenerations(obj.GenerationNumber + 1);
            end
            while ~isempty(child.AcceptedVoxelIndices)
                obj.AcceptedVoxelIndices{end + 1} = child.AcceptedVoxelIndices{1};
                child.AcceptedVoxelIndices(1) = [];
            end
            while ~isempty(child.RejectedVoxelIndices)
                obj.RejectedVoxelIndices{end + 1} = child.RejectedVoxelIndices{1};
                child.RejectedVoxelIndices(1) = [];
            end

        end

        function PruneAcceptedVoxels(obj)
            pruned_voxels = [];
            for index = 1 : length(obj.AcceptedVoxelIndices)
                next_voxel_set = obj.AcceptedVoxelIndices{index};
                size_voxel_set = numel(next_voxel_set);
                if (index == 1) || (size_voxel_set <= last_size)
                    pruned_voxels{index} = next_voxel_set;
                    last_size = size_voxel_set;
                else
                    obj.AcceptedVoxelIndices = pruned_voxels;
                    return;
                end
            end

        end

    end

    methods (Access = private)
        function AcceptAllPendingVoxelIndices(obj)
            while ~isempty(obj.PendingVoxelIndices)
                obj.AcceptedVoxelIndices{end + 1} = obj.PendingVoxelIndices{1};
                obj.PendingVoxelIndices(1) = [];
            end
        end

        function RejectAllPendingVoxelIndices(obj)
            while ~isempty(obj.PendingVoxelIndices)
                obj.RejectedVoxelIndices{end + 1} = obj.PendingVoxelIndices{1};
                obj.PendingVoxelIndices(1) = [];
            end
        end
    end

    methods (Static, Access = private)
        function concatenated_voxels = ConcatenateVoxels(voxels)
            concatenated_voxels = [];
            number_layers = length(voxels);
            for index = 1 : number_layers
                next_voxels = voxels{index};
                if ~isempty(next_voxels)
                    concatenated_voxels = cat(1, concatenated_voxels, next_voxels);
                end
            end
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function x = tinv(p,v);
% TINV   Inverse of Student's T cumulative distribution function (cdf).
%   X=TINV(P,V) returns the inverse of Student's T cdf with V degrees
%   of freedom, at the values in P.
%
%   The size of X is the common size of P and V. A scalar input
%   functions as a constant matrix of the same size as the other input.
%
% This is an open source function that was assembled by Eric Maris using
% open source subfunctions found on the web.

if nargin < 2,
    error('Requires two input arguments.');
end

[errorcode p v] = distchck(2,p,v);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

% Initialize X to zero.
x=zeros(size(p));

k = find(v < 0  | v ~= round(v));
if any(k)
    tmp  = NaN;
    x(k) = tmp(ones(size(k)));
end

k = find(v == 1);
if any(k)
  x(k) = tan(pi * (p(k) - 0.5));
end

% The inverse cdf of 0 is -Inf, and the inverse cdf of 1 is Inf.
k0 = find(p == 0);
if any(k0)
    tmp   = Inf;
    x(k0) = -tmp(ones(size(k0)));
end
k1 = find(p ==1);
if any(k1)
    tmp   = Inf;
    x(k1) = tmp(ones(size(k1)));
end

k = find(p >= 0.5 & p < 1);
if any(k)
    z = betainv(2*(1-p(k)),v(k)/2,0.5);
    x(k) = sqrt(v(k) ./ z - v(k));
end

k = find(p < 0.5 & p > 0);
if any(k)
    z = betainv(2*(p(k)),v(k)/2,0.5);
    x(k) = -sqrt(v(k) ./ z - v(k));
end

%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION distchck
%%%%%%%%%%%%%%%%%%%%%%%%%

function [errorcode,varargout] = distchck(nparms,varargin)
%DISTCHCK Checks the argument list for the probability functions.

errorcode = 0;
varargout = varargin;

if nparms == 1
    return;
end

% Get size of each input, check for scalars, copy to output
isscalar = (cellfun('prodofsize',varargin) == 1);

% Done if all inputs are scalars.  Otherwise fetch their common size.
if (all(isscalar)), return; end

n = nparms;

for j=1:n
   sz{j} = size(varargin{j});
end
t = sz(~isscalar);
size1 = t{1};

% Scalars receive this size.  Other arrays must have the proper size.
for j=1:n
   sizej = sz{j};
   if (isscalar(j))
      t = zeros(size1);
      t(:) = varargin{j};
      varargout{j} = t;
   elseif (~isequal(sizej,size1))
      errorcode = 1;
      return;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION betainv
%%%%%%%%%%%%%%%%%%%%%%%%%%%

function x = betainv(p,a,b);
%BETAINV Inverse of the beta cumulative distribution function (cdf).
%   X = BETAINV(P,A,B) returns the inverse of the beta cdf with
%   parameters A and B at the values in P.
%
%   The size of X is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   BETAINV uses Newton's method to converge to the solution.

%   Reference:
%      [1]     M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964.

%   B.A. Jones 1-12-93

if nargin < 3,
    error('Requires three input arguments.');
end

[errorcode p a b] = distchck(3,p,a,b);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

%   Initialize x to zero.
x = zeros(size(p));

%   Return NaN if the arguments are outside their respective limits.
k = find(p < 0 | p > 1 | a <= 0 | b <= 0);
if any(k),
   tmp = NaN;
   x(k) = tmp(ones(size(k)));
end

% The inverse cdf of 0 is 0, and the inverse cdf of 1 is 1.
k0 = find(p == 0 & a > 0 & b > 0);
if any(k0),
    x(k0) = zeros(size(k0));
end

k1 = find(p==1);
if any(k1),
    x(k1) = ones(size(k1));
end

% Newton's Method.
% Permit no more than count_limit interations.
count_limit = 100;
count = 0;

k = find(p > 0 & p < 1 & a > 0 & b > 0);
pk = p(k);

%   Use the mean as a starting guess.
xk = a(k) ./ (a(k) + b(k));

% Move starting values away from the boundaries.
if xk == 0,
    xk = sqrt(eps);
end
if xk == 1,
    xk = 1 - sqrt(eps);
end

h = ones(size(pk));
crit = sqrt(eps);

% Break out of the iteration loop for the following:
%  1) The last update is very small (compared to x).
%  2) The last update is very small (compared to 100*eps).
%  3) There are more than 100 iterations. This should NEVER happen.

while(any(abs(h) > crit * abs(xk)) & max(abs(h)) > crit    ...
                                 & count < count_limit),

    count = count+1;
    h = (betacdf(xk,a(k),b(k)) - pk) ./ betapdf(xk,a(k),b(k));
    xnew = xk - h;

% Make sure that the values stay inside the bounds.
% Initially, Newton's Method may take big steps.
    ksmall = find(xnew < 0);
    klarge = find(xnew > 1);
    if any(ksmall) | any(klarge)
        xnew(ksmall) = xk(ksmall) /10;
        xnew(klarge) = 1 - (1 - xk(klarge))/10;
    end

    xk = xnew;
end

% Return the converged value(s).
x(k) = xk;

if count==count_limit,
    fprintf('\nWarning: BETAINV did not converge.\n');
    str = 'The last step was:  ';
    outstr = sprintf([str,'%13.8f'],h);
    fprintf(outstr);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION betapdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = betapdf(x,a,b)
%BETAPDF Beta probability density function.
%   Y = BETAPDF(X,A,B) returns the beta probability density
%   function with parameters A and B at the values in X.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.

%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 26.1.33.

if nargin < 3,
   error('Requires three input arguments.');
end

[errorcode x a b] = distchck(3,x,a,b);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

% Initialize Y to zero.
y = zeros(size(x));

% Return NaN for parameter values outside their respective limits.
k1 = find(a <= 0 | b <= 0 | x < 0 | x > 1);
if any(k1)
   tmp = NaN;
    y(k1) = tmp(ones(size(k1)));
end

% Return Inf for x = 0 and a < 1 or x = 1 and b < 1.
% Required for non-IEEE machines.
k2 = find((x == 0 & a < 1) | (x == 1 & b < 1));
if any(k2)
   tmp = Inf;
    y(k2) = tmp(ones(size(k2)));
end

% Return the beta density function for valid parameters.
k = find(~(a <= 0 | b <= 0 | x <= 0 | x >= 1));
if any(k)
    y(k) = x(k) .^ (a(k) - 1) .* (1 - x(k)) .^ (b(k) - 1) ./ beta(a(k),b(k));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION betacdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function p = betacdf(x,a,b);
%BETACDF Beta cumulative distribution function.
%   P = BETACDF(X,A,B) returns the beta cumulative distribution
%   function with parameters A and B at the values in X.
%
%   The size of P is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   BETAINC does the computational work.

%   Reference:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 26.5.

if nargin < 3,
   error('Requires three input arguments.');
end

[errorcode x a b] = distchck(3,x,a,b);

if errorcode > 0
   error('Requires non-scalar arguments to match in size.');
end

% Initialize P to 0.
p = zeros(size(x));

k1 = find(a<=0 | b<=0);
if any(k1)
   tmp = NaN;
   p(k1) = tmp(ones(size(k1)));
end

% If is X >= 1 the cdf of X is 1.
k2 = find(x >= 1);
if any(k2)
   p(k2) = ones(size(k2));
end

k = find(x > 0 & x < 1 & a > 0 & b > 0);
if any(k)
   p(k) = betainc(x(k),a(k),b(k));
end

% Make sure that round-off errors never make P greater than 1.
k = find(p > 1);
p(k) = ones(size(k));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function options = bonminset(varargin)
%BONMINSET  Create or alter the options for Optimization with BONMIN
%
% options = bonminset('param1',value1,'param2',value2,...) creates an IPOPT
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the IPOPT
% default.
%
% options = bonminset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = bonminset() creates an options structure with all fields set to
% bonminset defaults.
%
% bonminset() prints a list of all possible fields and their function.
%
% See supplied BONMIN and IPOPT Documentation for further details of these options.

%   Copyright (C) 2011 Jonathan Currie (IPL)

% Print out possible values of properties.
if((nargin == 0) && (nargout == 0))
    printfields();
    return
end

%Names and Defaults
Names = {'algorithm','var_lin','cons_lin','allowable_fraction_gap','allowable_gap','cutoff','cutoff_decr','enable_dynamic_nlp','iteration_limit',...
         'nlp_failure_behavior','node_comparison','num_cut_passes','num_cut_passes_at_root',...
         'number_before_trust','number_strong_branch','solution_limit','sos_constraints','tree_search_strategy',...
         'variable_selection','feasibility_pump_objective_norm','heuristic_RINS','heuristic_dive_MIP_fractional',...
         'heuristic_dive_MIP_vectorLength','heuristic_dive_fractional','heuristic_dive_vectorLength',...
         'heuristic_feasibility_pump','pump_for_minlp','candidate_sort_criterion','maxmin_crit_have_sol',...
         'maxmin_crit_no_sol','min_number_strong_branch','number_before_trust_list','number_look_ahead',...
         'number_strong_branch_root','setup_pseudo_frac','coeff_var_threshold','dynamic_def_cutoff_decr','first_perc_for_cutoff_decr',...
         'max_consecutive_infeasible','num_resolve_at_infeasibles',...
         'num_resolve_at_node','num_resolve_at_root','second_perc_for_cutoff_decr','milp_solver','milp_strategy','bb_log_interval',...
         'bb_log_level','lp_log_level','milp_log_level','nlp_log_level','nlp_log_at_root','oa_log_frequency','oa_log_level',...
         'max_consecutive_failures','max_random_point_radius','num_retry_unsolved_random_point','random_point_perturbation_interval',...
         'random_point_type','warm_start'}';
Defaults = {'B-BB',[],[],0,0,1e100,1e-5,'no',2^31-1,'stop','best-bound',1,20,8,20,2^31-1,'enable','probed-dive','osi-strong',1,'no',...
            'yes','no','no','no','no','no','best-ps-cost',0.1,0.7,0,0,0,2^31-1,0.5,0.1,'no',-0.02,0,0,0,0,-0.05,'Cbc_D','find_good_sol',...
            100,1,0,0,0,0,100,1,10,100000,0,1,'Jon','none'}';
%Get IPOPT Defaults if required
if(nargin == 2 && strcmpi(varargin{2},'noipopt'))
    varargin = varargin(1);
else
    iset = ipoptset('bonmin');
    Names = [Names;fieldnames(iset)];
    Defaults = [Defaults;struct2cell(iset)];
end

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case {'first_perc_for_cutoff_decr','second_perc_for_cutoff_decr'}
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar non negative double
    case {'min_number_strong_branch','number_look_ahead','number_strong_branch_root','coeff_var_threshold'}
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar non zero double
    case {'oa_log_frequency','max_random_point_radius','random_point_perturbation_interval'}
        err = opticheckval.checkScalarGrtZ(value,field);
    %Scalar non negative integer
    case {'iteration_limit','num_cut_passes','num_cut_passes_at_root','number_before_trust','number_strong_branch',...
          'solution_limit','max_consecutive_infeasible','num_resolve_at_infeasibles','num_resolve_at_node',...
          'num_resolve_at_root','bb_log_interval','max_consecutive_failures','num_retry_unsolved_random_point'}
        err = opticheckval.checkScalarIntNonNeg(value,field);
    %Scalar double with bounds
    case 'cutoff_decr'
        err = opticheckval.checkScalarBoundLELE(value,field,-1e10,1e10);
    case {'allowable_fraction_gap','allowable_gap'}
        err = opticheckval.checkScalarBoundLELE(value,field,-1e20,1e20);
    case 'cutoff'
        err = opticheckval.checkScalarBoundLELE(value,field,-1e100,1e100);
    case {'maxmin_crit_have_sol','maxmin_crit_no_sol','setup_pseudo_frac'}
        err = opticheckval.checkScalarBoundLELE(value,field,0,1);
    case 'number_before_trust_list'
        err = opticheckval.checkScalarBoundLEL(value,field,-1,inf);
    %Scalar integer with bounds
    case 'feasibility_pump_objective_norm'
        err = opticheckval.checkScalarIntBoundLELE(value,field,1,2);
    case 'bb_log_level'
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,5);
    case {'lp_log_level','milp_log_level'}
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,4);
    case {'nlp_log_level','oa_log_level'}
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,2);
    case 'nlp_log_at_root'
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,12);
    %Vector integer 0/1
    case {'var_lin','cons_lin'}
        err = opticheckval.checkVector01(value,field);
    %Yes / No
    case {'heuristic_rins','heuristic_dive_mip_fractional','heuristic_dive_mip_vectorlength','heuristic_dive_fractional',...
          'heuristic_dive_vectorlength','heuristic_feasibility_pump','pump_for_minlp','dynamic_def_cutoff_decr'}
        err = opticheckval.checkYesNo(value,field);
    %Yes / No with Yes Crash
    case 'enable_dynamic_nlp'
        err = opticheckval.checkYesCrashNo(value,field);
    %Enable / Disable
    case 'sos_constraints'
        err = opticheckval.checkEnDis(value,field);
    %Misc String methods
    case 'algorithm'
        err = opticheckval.checkValidString(value,field,{'B-BB','B-OA','B-QG','B-Hyb','B-Ecp','B-iFP'});
    case 'nlp_failure_behavior'
        err = opticheckval.checkValidString(value,field,{'stop','fathom'});
    case 'node_comparison'
        err = opticheckval.checkValidString(value,field,{'best-bound','depth-first','breadth-first','dynamic','best-guess'});
    case 'tree_search_strategy'
        err = opticheckval.checkValidString(value,field,{'top-node','dive','probed-dive','dfs-dive','dfs-dive-dynamic'});
    case 'variable_selection'
        err = opticheckval.checkValidString(value,field,{'most-fractional','strong-branching','reliability-branching','qp-strong-branching',...
                                                         'lp-strong-branching','nlp-strong-branching','osi-simple','osi-strong','random'});
    case 'candidate_sort_criterion'
        err = opticheckval.checkValidString(value,field,{'best-ps-cost','worst-ps-cost','most-fractional','least-fractional'});
    case 'milp_solver'
        err = opticheckval.checkValidString(value,field,{'cbc_d','cplex'});
    case 'milp_strategy'
        err = opticheckval.checkValidString(value,field,{'find_good_sol','solve_to_optimality'});
    case 'random_point_type'
        err = opticheckval.checkValidString(value,field,{'jon','andreas','claudia'});
    case 'warm_start'
        err = opticheckval.checkValidString(value,field,{'none','fake_basis','optimum','interior_point'});

    %IPOPT or Unknown Parameter
    otherwise
        err = [];
        try
            %try for ipoptset arg
            iset = ipoptset(field,value);
        catch ME
            throw(ME);
        end
        if(isempty(iset))
            err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
        end
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf(' ----- BONMIN (MIP SOLVER) SETTINGS ----- \n');
fprintf('  BONMIN GENERAL SETTINGS:\n');
fprintf('                          algorithm: [ MIP Algorithm {''B-BB''}, ''B-OA'', ''B-QG'', ''B-Hyb'', ''B-Ecp'', ''B-iFP'' ] \n');
fprintf('                            var_lin: [ Decision Variable Linearity (0 - Nonlinear, 1 - Linear) {zeros(n,1)} ] \n');
fprintf('                           cons_lin: [ Constraint Linearity (0 - Nonlinear, 1 - Linear) {zeros(m,1)} ] \n');

fprintf('\n  BONMIN BRANCH & BOUND SETTINGS:\n');
fprintf('             allowable_fraction_gap: [ Value of Relative Gap when algorithm stops {0} ] \n');
fprintf('                      allowable_gap: [ Value of Absolute Gap when algorithm stops {0} ] \n');
fprintf('                             cutoff: [ The algorithm will only look for values better than cutoff {1e100} ] \n');
fprintf('                        cutoff_decr: [ Amount by which cutoff is decremented below a new best upper bound {1e-5} ] \n');
fprintf('                 enable_dynamic_nlp: [ Enable dynamic linear and quadratic rows addition in NLP ''yes'', {''no''} ] \n');
fprintf('                    iteration_limit: [ Cumulated maximum number of iterations in the algorithm to process continuous relaxations {2^31-1} ] \n');
fprintf('               nlp_failure_behavior: [ Action to take when an NLP is unsolved by IPOPT {''stop''}, ''fathom'' ] \n');
fprintf('                    node_comparison: [ Node selection strategy {''best-bound''}, ''depth-first'', ''breadth-first'', ''dynamic'', ''best-guess'' ] \n');
fprintf('                     num_cut_passes: [ Maximum number of cut passes at regular nodes {1} ] \n');
fprintf('             num_cut_passes_at_root: [ Maximum number of cut passes at root node {20} ] \n');
fprintf('                number_before_trust: [ Number of branches on a variable before its pseudo costs are to be believed in dynamic strong branching {8} ] \n');
fprintf('               number_strong_branch: [ Maximum number of variables considered for strong branching {20} ] \n');
fprintf('                     solution_limit: [ Abort after this many integer feasible solutions have been found (or 0 to deactivate) {2^31-1} ] \n');
fprintf('                    sos_constraints: [ Activate SOS Type 1 constraints {''enable''}, ''disable'' ] \n');
fprintf('               tree_search_strategy: [ Tree traversing strategy {''probed-dive''}, ''top-node'', ''dive'', ''dfs-dive'', ''dfs-dive-dynamic'' ] \n');
fprintf('                 variable_selection: [ Variable selection strategy {''osi-strong''}, ''most-fractional'', ''reliability-branching'', ''qp-strong-branching'', ''lp-strong-branching'', ''nlp-strong-branching'', ''osi-simple'', ''strong-branching'', ''random'' ] \n');

fprintf('\n  BONMIN MINLP HEURISTIC SETTINGS:\n');
fprintf('    feasibility_pump_objective_norm: [ Norm of feasibility pump objective norm {1}, 2 ] \n');
fprintf('                     heuristic_RINS: [ Use RINS heuristic {''no''}, ''yes'' ] \n');
fprintf('      heuristic_dive_MIP_fractional: [ Use Dive MIP fractional heuristic {''no''}, ''yes'' ] \n');
fprintf('    heuristic_dive_MIP_vectorLength: [ Use Dive MIP vectorLength heuristic {''no''}, ''yes'' ] \n');
fprintf('          heuristic_dive_fractional: [ Use Dive fractional heuristic ''no'', {''yes''} ] \n');
fprintf('        heuristic_dive_vectorLength: [ Use Dive vectorLength heuristic ''no'', {''yes''} ] \n');
fprintf('         heuristic_feasibility_pump: [ Use feasibility pump heuristic {''no''}, ''yes'' ] \n');
fprintf('                     pump_for_minlp: [ Use FP for MINLP {''no''}, ''yes'' ] \n');

fprintf('\n  BONMIN NON-CONVEX PROBLEM SETTINGS:\n');
fprintf('                coeff_var_threshold: [ Coefficient of variation threshold (for dynamic definition of cutoff_decr) {0.1} ] \n');
fprintf('            dynamic_def_cutoff_decr: [ Define cutoff_decr dynamically? {''no''}, ''yes'' ] \n');
fprintf('         first_perc_for_cutoff_decr: [ Percentage used when the coeff of variance is smaller than the threshold, to compute cutoff_decr dynamically {-0.02} ] \n');
fprintf('         max_consecutive_infeasible: [ Maximum number of consecutive infeasible subproblems before aborting a branch {0} ] \n');
fprintf('         num_resolve_at_infeasibles: [ Number of tries to resolve an infeasible node (other than the root) of the tree with different starting points {0} ] \n');
fprintf('                num_resolve_at_node: [ Number of tries to resolve a node (other than the root) of the tree with different starting points (like multi-start) {0} ] \n');
fprintf('                num_resolve_at_root: [ Number of tries to resolve the root node with different starting points {0} ] \n');
fprintf('        second_perc_for_cutoff_decr: [ Percentage used when the coeff of variance is greater than the threshold, to compute cutoff_decr dynamically {-0.05} ] \n');

fprintf('\n  BONMIN NLP SOLUTION ROBUSTNESS SETTINGS:\n');
fprintf('           max_consecutive_failures: [ Number of consecutive unsolved problems before aborting a branch of the tree {10} ] \n');
fprintf('            max_random_point_radius: [ Max value for coordinate of a random point {100000} ] \n');
fprintf('    num_retry_unsolved_random_point: [ Number of times that the algorithm will try to resolve an unsolved NLP with a random starting point (with new point) {0} ] \n');
fprintf(' random_point_perturbation_interval: [ Amount by which starting point is perturbed when choosing a random point {1} ] \n');
fprintf('                  random_point_type: [ Random starting point method {''Jon''}, ''Andreas'', ''Claudia'' ] \n');
fprintf('                         warm_start: [ Warm start method {''none''}, ''fake_basis'', ''optimum'', ''interior'' ] \n');

fprintf('\n  BONMIN STRONG BRANCHING SETTINGS:\n');
fprintf('           candidate_sort_criterion: [ Criterion used to choose candidates in strong-branching {''best-ps-cost''}, ''worst-ps-cost'', ''most-fractional'', ''least-fractional'' ] \n')
fprintf('               maxmin_crit_have_sol: [ Weight towards minimum of lower and upper branching estimates when a solution has been found {0.1} ] \n')
fprintf('                 maxmin_crit_no_sol: [ Weight towards minimum of lower and upper branching estimates when no solution has been found {0.7} ] \n')
fprintf('           min_number_strong_branch: [ Minimum number of variables for strong branching (overriding trust) {0} ] \n')
fprintf('           number_before_trust_list: [ Number of branches on a variable vbefore its pseudo costs are to be believed during setup of strong branching list {0} ] \n')
fprintf('                  number_look_ahead: [ Limit of look-ahead strong-branching trials {0} ] \n')
fprintf('          number_strong_branch_root: [ Maximum number of variables considered for strong branching at root node {2^31-1} ] \n')
fprintf('                  setup_pseudo_frac: [ Proportion of strong branching list that has to be taken from most integer infeasible list {0.5} ] \n')

fprintf('\n  BONMIN MILP SOLVER SETTINGS:\n');
fprintf('                        milp_solver: [ Subsolver to solve MILP problems in OA decompositions {''Cbc_D''}, ''Cplex'' (must be installed on your PC) ] \n');
fprintf('                      milp_strategy: [ MILP solving strategy {''find_good_sol''}, ''solve_to_optimality'' ] \n');

fprintf('\n  BONMIN OUPUT AND LOG LEVEL SETTINGS:\n');
fprintf('                    bb_log_interval: [ Interval at which node level output is printed (number of nodes) {100} ] \n');
fprintf('                       bb_log_level: [ Level of branch and bound log detail (0-5) {1} ] \n');
fprintf('                       lp_log_level: [ Level of LP solver log detail (0-4) {0} ] \n');
fprintf('                     milp_log_level: [ Level of MILP solver log detail (0-4) {0} ] \n');
fprintf('                      nlp_log_level: [ Level of NLP solver log detail (0-2) {0} ] \n');
fprintf('                    nlp_log_at_root: [ Level of NLP solver log detail at root node (0-12) {0} ] \n');
fprintf('                   oa_log_frequency: [ Frequency (in seconds) of OA log messages  {100} ] \n');
fprintf('                       oa_log_level: [ Level of OA decomposition log detail (0-2) {1} ] \n');

fprintf('\n ----- IPOPT (RELAXED SOLVER) SETTINGS ----- \n');
ipoptset('bonmin');

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%#eml
function [min,med,max] = get_median_2d(inbuf)

[nrows ncols] = size(inbuf);

tbuf = inbuf;

for ii=eml.unroll(1:nrows)
    tbuf(ii, :) = get_median_1d(inbuf(ii, :));
end

for jj=eml.unroll(1:ncols)
    tbuf(:, jj) = get_median_1d(tbuf(:, jj)')';
end

min = tbuf(1, 1);
med = tbuf(ceil(nrows/2), ceil(ncols/2));
max = tbuf(nrows, ncols);

end

author: Sable

repo: mcbench-benchmarks


text:
begin=clock;
omega=1;
over_r=2.4; %note that over_r should be between 3 and 4 for this program.
rate=5.3; %this is approximatly the ratio N/r, here the ratio must be above 1
n=ceil(rate*over_r);%this is the number of unions of sets
p=ceil(2*n-over_r); %this is the degree of the fine mesh, use p=n+1, or larger
L=200/n; %the signal is sampled from -L*T:L*T
%L=0.5;
step=over_r/2/omega;

%R1=round(-2/3+n/3);
%L1=R1+1-n;
%L2=round(2/3-n/3);
%R2=n+L2-1;

%I have doubled over_r and halved L so that the
%interval is the same, i.e., width is the same
%as it was for the single sample case.

rand('state',sum(100*clock))
%rand('seed',4)
h=(rand(n,1)-1/2)*step;
%h=(-1/2:1/n:1/2-1/n)*step;

width=L*step;
y=-width:step/p:width;
ysize=size(y);
ysize=ysize(2);
% y is the sampling set for signal_origin and y+h for signal_shift
signal=zeros(n,ysize);
pad=zeros(n,ysize);

signal_origin=zeros(size(y));
signal_shift=zeros(size(y));
pad_origin=zeros(size(y));
pad_shift=zeros(size(y));

%rand('seed',0)
%this sets the seed to a fixed value, 0, so that I get reproducable results
rand('state',sum(100*clock))
R=100;
%the first two collumns are the amplitudes, real(first) and complex(second),
%i.e., signal_coef(1,q)+i*signal_coef(2,q), the third and fourth are the
%band limits, for example, min(signal_coef(3,q),signal_coef(4,q)) is
%the left bandwidth, the max will give the right bandwidth
signal_coef=rand(4,R);
signal_coef(1,:)=2*(signal_coef(1,:)-1/2);
signal_coef(1,:)=signal_coef(1,:)/norm(signal_coef(1,:),2)/2/pi;
signal_coef(2,:)=2*(signal_coef(2,:)-1/2);
signal_coef(2,:)=signal_coef(2,:)/norm(signal_coef(2,:),2)/2/pi;
signal_coef(3,:)=2*(signal_coef(3,:)-1/2);
signal_coef(3,:)=omega*signal_coef(3,:)/max(abs(signal_coef(3,:)));
signal_coef(4,:)=2*(signal_coef(4,:)-1/2);
signal_coef(4,:)=omega*signal_coef(4,:)/max(abs(signal_coef(4,:)));

%need an axis for the dual space that has the same number of
%elements as y does.

dual_axis=2*p*omega/over_r/(2*p*L+1)*(-p*L:1:p*L);

sample_dual=zeros(size(dual_axis));
dualsize=size(dual_axis);
dualsize=dualsize(2);

for q=1:R
  for j=1:dualsize
if dual_axis(j)>=min(signal_coef(3,q),signal_coef(4,q)) & dual_axis(j)<=max(signal_coef(3,q),signal_coef(4,q))
      sample_dual(j)=sample_dual(j)+signal_coef(1,q)+i*signal_coef(2,q);
    end
  end
end

tmp=0;
for q=1:R
  signal_origin=signal_origin+(signal_coef(1,q)+i*signal_coef(2,q))./(sqrt(2*pi)*2*pi*i*y).*(exp(2*pi*i*y*max(signal_coef(3,q),signal_coef(4,q)))-exp(2*pi*i*y*min(signal_coef(3,q),signal_coef(4,q))))*2*pi;
%  signal_shift=signal_shift+(signal_coef(1,q)+i*signal_coef(2,q))./(sqrt(2*pi)*2*pi*i*(y+h(2))).*(exp(2*pi*i*(y+h(2))*max(signal_coef(3,q),signal_coef(4,q)))-exp(2*pi*i*(y+h(2))*min(signal_coef(3,q),signal_coef(4,q))))*2*pi;
  tmp=tmp+(signal_coef(1,q)+i*signal_coef(2,q))*(max(signal_coef(3,q),signal_coef(4,q))-min(signal_coef(3,q),signal_coef(4,q)))/sqrt(2*pi)*2*pi;
end
  joe=0;
  for q=1:ysize
    if y(q)==0
      joe=q;
    end
  end
  if joe>0
    signal_origin(joe)=tmp;
  end

for j=1:n
tmp=0;
for q=1:R
  signal(j,:)=signal(j,:)+(signal_coef(1,q)+i*signal_coef(2,q))./(sqrt(2*pi)*2*pi*i*(y+h(j))).*(exp(2*pi*i*(y+h(j))*max(signal_coef(3,q),signal_coef(4,q)))-exp(2*pi*i*(y+h(j))*min(signal_coef(3,q),signal_coef(4,q))))*2*pi;
  tmp=tmp+(signal_coef(1,q)+i*signal_coef(2,q))*(max(signal_coef(3,q),signal_coef(4,q))-min(signal_coef(3,q),signal_coef(4,q)))/sqrt(2*pi)*2*pi;
end
  joe=0;
  for q=1:ysize
    if y(q)+h(j)==0
      joe=q;
    end
  end
  if joe>0
    signal(j,joe)=tmp;
  end
end

% p=n is the measure of difference between the fine and coarse mesh
for j=1:p:ysize
  pad(:,j)=signal(:,j);
end

pad_dual=zeros(n,ysize);

for j=1:dualsize
  for k=1:n
    pad_dual(k,j)=sum(pad(k,:).*exp(-2*pi*i*dual_axis(j)*y));
  end
end

pad_dual=pad_dual*(over_r/2/omega/p)/sqrt(2*pi);

%SHOULD THEY BE MULTIPLIED BY p SO THAT THEY MATCH sample_dual?

%this makes it so that the origin matches.
for j=1:n
  pad_dual(j,:)=pad_dual(j,:).*exp(-2*pi*i*h(j)*dual_axis);
end

kappa=max(2,min(n,floor((n+over_r+1)/(n-over_r+1))))

%b=zeros(2*kappa,1);
%B=zeros(2*kappa);
%B(1,1)=1; B(1,2)=1; B(1,3)=-1;
%B(kappa,end)=1; B(kappa,end-1)=1; B(kappa,end-2)=-1;
%
%for j=2:kappa-1
%  B(j,2*(j-1))=-1;
%  B(j,2*(j-1)+1)=1;
%  B(j,2*(j-1)+2)=1;
%  B(j,2*(j-1)+3)=-1;
%end
%b(1)=-1;
%b(kappa)=1;
%for j=1:kappa
%  B(j+kappa,2*j-1)=-1;
%  B(j+kappa,2*j)=1;
%  b(j+kappa)=n-1;
%end
%tmp1=B\b;

zones=zeros(kappa,2);
for j=1:kappa
%zones(j,1)=tmp1(2*j-1);
%zones(j,2)=tmp1(2*j);
zones(j,1)=round(j*(n+1)/(kappa+1)-n);
zones(j,2)=zones(j,1)+n-1;
end

%for j=1:kappa
%  if abs(zones(j,1)-round(zones(j,1)))<=abs(zones(j,2)-round(zones(j,2)))
%    zones(j,1)=round(zones(j,1));
%    zones(j,2)=n-1+zones(j,1);
%  else
%    zones(j,2)=round(zones(j,2));
%    zones(j,1)=zones(j,2)-n+1;
%  end
%end

A=zeros(n);
for j=1:n
  for k=1:n
    A(j,k)=exp(2*pi*i*h(k)*j/step);
  end
end
A_inverse=inv(A);

c=zeros(kappa,n);

for j=1:kappa
  e=zeros(n,1);
  e(1-zones(j,1))=1;
  R=diag(exp(2*pi*i*h*(zones(j,1)-1)/step));
  c(j,:)=(inv(R)*A_inverse*e)'; %'
end

%e=zeros(n,1);
%e(1-L2)=1;
%R=diag(exp(2*pi*i*h*(L2-1)/step));
%c(2,:)=(inv(R)*A_inverse*e)'; %'

partition=over_partitions(zones,omega,over_r,n,dual_axis);

%left_partition=zeros(size(dual_axis));
%right_partition=zeros(size(dual_axis));

%for j=1:length(dual_axis)
%  if dual_axis(j)<=(L1-1)/step+omega
%    left_partition(j)=0;
%  elseif (L1-1)/step+omega<dual_axis(j) & dual_axis(j)<-omega
%    left_partition(j)=rho((-dual_axis(j)-omega)/(-(L1-1)/step-2*omega));
%  elseif -omega<=dual_axis(j) & dual_axis(j)<=(L2-1)/step+omega
%    left_partition(j)=1;
%  elseif (L2-1)/step+omega<dual_axis(j) & dual_axis(j)<(R1+1)/step-omega
%    left_partition(j)=rho((dual_axis(j)-((L2-1)/step+omega))/((R1-L2+2)/step-2*omega));
%  else
%    left_partition(j)=0;
%  end
%end
%for j=1:length(dual_axis)
%  if dual_axis(j)<=(L2-1)/step+omega
%    right_partition(j)=0;
%  elseif (L2-1)/step+omega<dual_axis(j) & dual_axis(j)<(R1+1)/step-omega
%    right_partition(j)=1-rho((dual_axis(j)-((L2-1)/step+omega))/((R1-L2+2)/step-2*omega));
%  elseif (R1+1)/step-omega<=dual_axis(j) & dual_axis(j)<=omega
%    right_partition(j)=1;
%  elseif omega<dual_axis(j) & dual_axis(j)<(R2+1)/step-omega
%    right_partition(j)=rho((dual_axis(j)-omega)/((R2+1)/step-2*omega));
%  else
%    right_partition(j)=0;
%  end
%end

filters=zeros(n,length(dual_axis));
for j=1:n
  for k=1:kappa
    filters(j,:)=filters(j,:)+c(k,j)*partition(k,:);
  end
end

filtered_duals=zeros(size(filters));

for j=1:n
  filtered_duals(j,:)=filtered_duals(j,:)+filters(j,:).*pad_dual(j,:);
end

recon_dual=zeros(size(dual_axis));
for j=1:n
  recon_dual=recon_dual+filtered_duals(j,:);
end

recon_signal=zeros(size(y));

k=-p*L:1:p*L;

for j=1:dualsize
  recon_signal=recon_signal+recon_dual(j)*exp(2*pi*i*k/(2*p*L+1)*k(j));
end

recon_signal=recon_signal/(2*p*L+1)*sqrt(2*pi)*omega*p/over_r*2*p;

%hold off
semilogy(y,abs(recon_signal-signal_origin),'k');
%pause
%hold on
%plot(h,10^(-2),'o');

%h/step
cond(A)

done=clock;
time_cost=done-begin;
time_cost=time_cost(6)+60*time_cost(5)+360*time_cost(4)

%L
%step

%for j=1:n
%  hold off
%  plot(dual_axis,real(filters(j,:)))
%  hold on
%  plot(dual_axis,imag(filters(j,:)),'r')
%  hold off
%  pause
%end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [ mX ] = ProjectSymmetricMatrixSet( mX )
% ----------------------------------------------------------------------------------------------- %
% [ mX ] = ProjectSymmetricMatrixSet( mX )
%   Projecting the input matrix into the Convex Set of Symmetric Matrices.
% Input:
%   - mX            -   Input Matrix.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - mX            -   Output Matrix.
%                       Symmteirc Matrix which is the orthogonal projection
%                       of the input matrix.
%                       Structure: Matrix.
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  h
% Remarks:
%   1.  a
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     15/08/2018  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

mX = (mX.' + mX) / 2;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
% this script tests 1st Order Statistics features between CERR and pyradiomics.
%
% RKP, 03/22/2018

firstOrderParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','tests_for_cerr','test_first_order_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','data_for_cerr_tests','CERR_plans','head_neck_ex1_20may03.mat.bz2');

planC = loadPlanC(cerrFileName,tempdir);
indexS = planC{end};

paramS = getRadiomicsParamTemplate(firstOrderParamFileName);
strNum = getMatchingIndex(paramS.structuresC{1},{planC{indexS.structures}.structureName});
scanNum = getStructureAssociatedScan(strNum,planC);

%% Calculate features using CERR

firstOrderS = calcGlobalRadiomicsFeatures...
            (scanNum, strNum, paramS, planC);
firstOrderS = firstOrderS.Original.firstOrderS;
cerrFirstOrderV = [firstOrderS.energy, firstOrderS.totalEnergy, firstOrderS.interQuartileRange, ...
    firstOrderS.kurtosis+3, firstOrderS.max, firstOrderS.mean, firstOrderS.meanAbsDev, ...
    firstOrderS.median, firstOrderS.medianAbsDev, firstOrderS.min, ...
    firstOrderS.P10, firstOrderS.P90, firstOrderS.interQuartileRange, ...
    firstOrderS.robustMeanAbsDev, firstOrderS.rms, firstOrderS.skewness, ...
    firstOrderS.std, firstOrderS.var, firstOrderS.entropy];

%% Calculate features using pyradiomics

testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
    single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
mask3M = zeros(size(testM),'logical');
[rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
[maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
mask3M(:,:,uniqueSlices) = maskBoundBox3M;

scanType = 'original';

dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
pixelSize = [dx dy dz]*10;

teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);

pyradFirstorderNamC = {'Energy', 'TotalEnergy','InterquartileRange','Kurtosis',...
    'Maximum', 'Mean','MeanAbsoluteDeviation','Median','medianAbsDev',...
    'Minimum','10Percentile','90Percentile','InterquartileRange',...
    'RobustMeanAbsoluteDeviation','RootMeanSquared','Skewness',...
    'StandardDeviation','Variance','Entropy'};

pyradFirstorderNamC = strcat(['original', '_firstorder_'],pyradFirstorderNamC);

pyRadFirstOrderV = [];
for i = 1:length(pyradFirstorderNamC)
    if isfield(teststruct,pyradFirstorderNamC{i})
        pyRadFirstOrderV(i) = teststruct.(pyradFirstorderNamC{i});
    else
        pyRadFirstOrderV(i) = NaN;
    end
end

diffFirstOrderV = (cerrFirstOrderV - pyRadFirstOrderV) ./ cerrFirstOrderV * 100

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [R,W] = rodrigues(N,phi)
  % RODRIGUES Construct the rotation matrix for rotating vectors around each
  % given vector.
  %
  % Inputs:
  %   N  #N by 3 list of 3D vectors
  %   phi  angle to rotate by.
  % Outputs:
  %   R  3#N by 3#N "rotation" matrix
  %   W  3#N by 3#N cross product matrix
  %
  % Example:
  %   % Given a scalar function Z on a mesh (V,F)
  %   G = grad(V,F);
  %   N = normalizerow(normals(V,F));
  %   X = reshape(G*Z,[],3);
  %   Y = reshape(rodrigues(N,pi/2)*G*Z,[],3);
  %   BC = barycenter(V,F);
  %   tsurf(F,V,'CData',Z,fphong,'EdgeColor','none');
  %   hold on;
  %   quiver3(BC(:,1),BC(:,2),BC(:,3),X(:,1),X(:,2),X(:,3))
  %   quiver3(BC(:,1),BC(:,2),BC(:,3),Y(:,1),Y(:,2),Y(:,3))
  %   hold off;
  %

  m = size(N,1);
  % Either uniform angle or per-vector
  if numel(phi)>1
    assert(numel(phi) == m);
    phi = diag(sparse(phi));
  end
  % https://math.stackexchange.com/a/142831
  %
  W = sparse( ...
    [1 2 0 2 0 1]*m+(1:m)', ...
    [0 0 1 1 2 2]*m+(1:m)', ...
    [N(:,3) -N(:,2) -N(:,3) N(:,1) N(:,2) -N(:,1)], ...
    3*m,3*m);
  if isempty(phi)
    R = [];
  else
    R = speye(3*m,3*m) + sin(phi)*W + (2*sin(phi/2)^2)*W*W;
  end
end

author: alecjacobson

repo: gptoolbox


text:
function divdif_test03 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST03 tests DIF_BASIS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 June 2013
%
%  Author:
%
%    John Burkardt
%
  ntab = 5;
  nstep = 9;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST03\n' );
  fprintf ( 1, '  DIF_BASIS computes Lagrange basis polynomials\n' );
  fprintf ( 1, '  in difference form.\n' );
  fprintf ( 1, '\n' );
%
%  Set the base points.
%
  xtab = r8vec_indicator ( ntab );

  r8vec_print ( ntab, xtab, '  The base points:' );
%
%  Get the difference tables for the basis polynomials and print them.
%
  diftab = dif_basis ( ntab, xtab );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The table of difference vectors defining the basis\n' );
  fprintf ( 1, '  polynomials.  Each column represents a polynomial.\n' );
  fprintf ( 1, '\n' );
  for i = 1 : ntab
    fprintf ( 1, '  ' );
    for j = 1 : ntab
      fprintf ( 1, '%14f', diftab(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Evaluate basis polynomial 3 at a set of points.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Evaluate basis polynomial #3 at a set of points.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X        Y\n' );
  fprintf ( 1, '\n' );
  xhi = ntab;
  xlo = 1.0;

  for i = 1 : nstep

    xval = ( ( nstep - i     ) * xlo   ...
           + (         i - 1 ) * xhi ) ...
           / ( nstep     - 1 );

    yval = dif_val ( ntab, xtab, diftab(1:ntab,3), xval );

    fprintf ( 1, '  %14f  %14f\n', xval, yval );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function vl_demo_alldist
%

numRepetitions  = 3 ;
numDimensions   = 1000 ;
numSamplesRange = [300] ;

settingsRange = {{'alldist2', 'double',  'l2',         }, ...
                 {'alldist',  'double',  'l2', 'nosimd'}, ...
                 {'alldist',  'double',  'l2'          }, ...
                 {'alldist2', 'single',  'l2',         }, ...
                 {'alldist',  'single',  'l2', 'nosimd'}, ...
                 {'alldist',  'single',  'l2'          }, ...
                 {'alldist2', 'double',  'l1',         }, ...
                 {'alldist',  'double',  'l1', 'nosimd'}, ...
                 {'alldist',  'double',  'l1'          }, ...
                 {'alldist2', 'single',  'l1',         }, ...
                 {'alldist',  'single',  'l1', 'nosimd'}, ...
                 {'alldist',  'single',  'l1'          }, ...
                 {'alldist2', 'double',  'chi2',         }, ...
                 {'alldist',  'double',  'chi2', 'nosimd'}, ...
                 {'alldist',  'double',  'chi2'          }, ...
                 {'alldist2', 'single',  'chi2',         }, ...
                 {'alldist',  'single',  'chi2', 'nosimd'}, ...
                 {'alldist',  'single',  'chi2'          }, ...
                 {'alldist2', 'double',  'hell',         }, ...
                 {'alldist',  'double',  'hell', 'nosimd'}, ...
                 {'alldist',  'double',  'hell'          }, ...
                 {'alldist2', 'single',  'hell',         }, ...
                 {'alldist',  'single',  'hell', 'nosimd'}, ...
                 {'alldist',  'single',  'hell'          }, ...
                 {'alldist2', 'double',  'kl2',         }, ...
                 {'alldist',  'double',  'kl2', 'nosimd'}, ...
                 {'alldist',  'double',  'kl2'          }, ...
                 {'alldist2', 'single',  'kl2',         }, ...
                 {'alldist',  'single',  'kl2', 'nosimd'}, ...
                 {'alldist',  'single',  'kl2'          }, ...
                 {'alldist2', 'double',  'kl1',         }, ...
                 {'alldist',  'double',  'kl1', 'nosimd'}, ...
                 {'alldist',  'double',  'kl1'          }, ...
                 {'alldist2', 'single',  'kl1',         }, ...
                 {'alldist',  'single',  'kl1', 'nosimd'}, ...
                 {'alldist',  'single',  'kl1'          }, ...
                 {'alldist2', 'double',  'kchi2',         }, ...
                 {'alldist',  'double',  'kchi2', 'nosimd'}, ...
                 {'alldist',  'double',  'kchi2'          }, ...
                 {'alldist2', 'single',  'kchi2',         }, ...
                 {'alldist',  'single',  'kchi2', 'nosimd'}, ...
                 {'alldist',  'single',  'kchi2'          }, ...
                 {'alldist2', 'double',  'khell',         }, ...
                 {'alldist',  'double',  'khell', 'nosimd'}, ...
                 {'alldist',  'double',  'khell'          }, ...
                 {'alldist2', 'single',  'khell',         }, ...
                 {'alldist',  'single',  'khell', 'nosimd'}, ...
                 {'alldist',  'single',  'khell'          }, ...
                } ;

%settingsRange = settingsRange(end-5:end) ;

styles = {} ;
for marker={'x','+','.','*','o'}
  for color={'r','g','b','k','y'}
    styles{end+1} = {'color', char(color), 'marker', char(marker)} ;
  end
end

for ni=1:length(numSamplesRange)
  for ti=1:length(settingsRange)
    tocs = [] ;
    for ri=1:numRepetitions
      rand('state',ri) ;
      randn('state',ri) ;
      numSamples = numSamplesRange(ni) ;
      settings = settingsRange{ti} ;
      [tocs(end+1), D] = run_experiment(numDimensions, ...
                                        numSamples, ...
                                        settings) ;
    end
    means(ni,ti) = mean(tocs) ;
    stds(ni,ti)  = std(tocs) ;
    if mod(ti-1,3) == 0
      D0 = D ;
    else
      err = max(abs(D(:)-D0(:))) ;
      fprintf('err %f\n', err) ;
      if err > 1, keyboard ; end
    end
  end
end

if 0
  figure(1) ; clf ; hold on ;
  numStyles = length(styles) ;
  for ti=1:length(settingsRange)
    si = mod(ti - 1, numStyles) + 1 ;
    h(ti) = plot(numSamplesRange, means(:,ti), styles{si}{:}) ;
    leg{ti} = sprintf('%s ', settingsRange{ti}{:}) ;
    errorbar(numSamplesRange, means(:,ti), stds(:,ti), 'linestyle', 'none') ;
  end
end

for ti=1:length(settingsRange)
  leg{ti} = sprintf('%s ', settingsRange{ti}{:}) ;
end

figure(1) ; clf ;
barh(means(end,:)) ;
set(gca,'ytick', 1:length(leg), 'yticklabel', leg,'ydir','reverse') ;
xlabel('Time [s]') ;

function [elaps, D] = run_experiment(numDimensions, numSamples, settings)

distType  = 'l2' ;
algType   = 'alldist' ;
classType = 'double' ;
useSimd   = true ;

for si=1:length(settings)
  arg = settings{si} ;
  switch arg
    case {'l1', 'l2', 'chi2', 'hell', 'kl2', 'kl1', 'kchi2', 'khell'}
      distType = arg ;
    case {'alldist', 'alldist2'}
      algType = arg ;
    case {'single', 'double'}
      classType = arg ;
    case 'simd'
      useSimd = true ;
    case 'nosimd'
      useSimd = false ;
    otherwise
      assert(false) ;
  end
end

X = rand(numDimensions, numSamples) ;
X(X < .3) = 0 ;

switch classType
  case 'double'
  case 'single'
    X = single(X) ;
end

vl_simdctrl(double(useSimd)) ;

switch algType
  case 'alldist'
    tic ; D = vl_alldist(X, distType) ; elaps = toc ;
  case 'alldist2'
    tic ; D = vl_alldist2(X, distType) ; elaps = toc ;
end

author: Cloud-CV

repo: object-proposals


text:
function add_header_to_bin(bin_filename, fc_requested, fc_programmed, fs_requested, fs_programmed)

[fc_requested_exist, fc_programmed_exist, fs_requested_exist, fs_programmed_exist] = read_header_from_bin(bin_filename);
if fc_requested_exist ~= inf
    disp('There is already a header! Just return.');
    disp(['fc_requested ' num2str(fc_requested_exist) ' fc_programmed ' num2str(fc_programmed_exist)  ' fs_requested ' num2str(fs_requested_exist)  ' fs_programmed ' num2str(fs_programmed_exist) ]);
    return;
end

fid = fopen(bin_filename, 'r');

if fid==-1
    disp('add_header_to_bin: Can not open file for read!');
    return;
end

tmp_store = fread(fid, inf, 'uint8');
fclose(fid);

fc_requested_magic = 73492.215;
fc_programmed_magic = -0.7923597;
fs_requested_magic = -189978508;
fs_programmed_magic = 93.126712;

reserve1_magic = -53243.129;
reserve2_magic = 0.0008123898;
reserve3_magic = -6.0098321;
reserve4_magic = 237.09983;

fid = fopen(bin_filename, 'w');

if fid==-1
    disp('add_header_to_bin: Can not open file for write!');
    return;
end

fwrite(fid, fc_requested_magic, 'double');
fwrite(fid, fc_requested, 'uint64');

fwrite(fid, fc_programmed_magic, 'double');
fwrite(fid, fc_programmed, 'uint64');

fwrite(fid, fs_requested_magic, 'double');
fwrite(fid, fs_requested, 'uint64');

fwrite(fid, fs_programmed_magic, 'double');
fwrite(fid, fs_programmed, 'uint64');

fwrite(fid, reserve1_magic, 'double');
fwrite(fid, 0, 'uint64');

fwrite(fid, reserve2_magic, 'double');
fwrite(fid, 0, 'uint64');

fwrite(fid, reserve3_magic, 'double');
fwrite(fid, 0, 'uint64');

fwrite(fid, reserve4_magic, 'double');
fwrite(fid, 0, 'uint64');

fwrite(fid, tmp_store, 'uint8');

fclose(fid);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function varargout = encriptacionprograma(varargin)
%ENCRIPTACIONPROGRAMA M-file for encriptacionprograma.fig
%      ENCRIPTACIONPROGRAMA, by itself, creates a new ENCRIPTACIONPROGRAMA or raises the existing
%      singleton*.
%
%      H = ENCRIPTACIONPROGRAMA returns the handle to a new ENCRIPTACIONPROGRAMA or the handle to
%      the existing singleton*.
%
%      ENCRIPTACIONPROGRAMA('Property','Value',...) creates a new ENCRIPTACIONPROGRAMA using the
%      given property value pairs. Unrecognized properties are passed via
%      varargin to encriptacionprograma_OpeningFcn.  This calling syntax produces a
%      warning when there is an existing singleton*.
%
%      ENCRIPTACIONPROGRAMA('CALLBACK') and ENCRIPTACIONPROGRAMA('CALLBACK',hObject,...) call the
%      local function named CALLBACK in ENCRIPTACIONPROGRAMA.M with the given input
%      arguments.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help encriptacionprograma

% Last Modified by GUIDE v2.5 11-Jun-2013 18:21:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @encriptacionprograma_OpeningFcn, ...
                   'gui_OutputFcn',  @encriptacionprograma_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 encriptacionprograma is made visible.
function encriptacionprograma_OpeningFcn(hObject, ~, 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   unrecognized PropertyName/PropertyValue pairs from the
%            command line (see VARARGIN)

% Choose default command line output for encriptacionprograma
handles.output = hObject;
axes(handles.axes5);
BackGr = imread('encrip.gif');
  imshow(BackGr);

% Update handles structure
guidata(hObject, handles);
clc;
global Img;
global EncImg;
global DecImg;
% UIWAIT makes encriptacionprograma wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = encriptacionprograma_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
varargout{1} = handles.output;
% --- Executes during object creation, after setting all properties.
function imagenin_CreateFcn(hObject, eventdata, handles)
% hObject    handle to imagenin (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate imagenin

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global Img;
global key;
X = uigetfile('*.jpg;*.tiff;*.ppm;*.pgm;*.png','pick a jpge file');
Img = imread(X);
  axes(handles.imagenin);
  imshow(Img);

[n m k] = size(Img);
  key = keyGen(n*m);
guidata(hObject, handles);
function [proImageOut] = imageProcess(ImgInp,key)
[n m k] = size(ImgInp);
% key =cell2mat(struct2cell( load('key5.mat')));
% key = keyGen(n*m);
for ind = 1 : m
    Fkey(:,ind) = key((1+(ind-1)*n) : (((ind-1)*n)+n));
end
len = n;
bre = m;
for ind = 1 : k
    Img = ImgInp(:,:,ind);
for ind1 = 1 : len
    for ind2 = 1 : bre
        proImage(ind1,ind2) = bitxor(Img(ind1,ind2),Fkey(ind1,ind2));
    end
end
proImageOut(:,:,ind) = proImage(:,:,1);
end
figure,imshow(proImageOut);
return;

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global Img ;
global EncImg;
global key;
EncImg = imageProcess(Img,key);
axes(handles.axes2)
imshow(EncImg);
imwrite(EncImg,'nueva imagen encriptada.jpg','jpg');

guidata(hObject, handles);

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global DecImg;
global EncImg;
global key;
DecImg = imageProcess(EncImg,key);
axes(handles.axes3);
imshow(DecImg);
imwrite(DecImg,'imagen decriptada.jpg','jpg');
guidata(hObject, handles);
function [key] = keyGen(n)
n = n*4;
% n = 2048*2048*16;
%n = 24 * 24 * 8;
bin_x = zeros(n,1,'uint8');
r = 3.9999998;
bin_x_N_Minus_1 =  0.300001;
x_N = 0;
tic
for ind = 2 : n
    x_N = 1 - 2* bin_x_N_Minus_1 * bin_x_N_Minus_1;
     if (x_N > 0.0)
        bin_x(ind-1) = 1;
    end
     bin_x_N_Minus_1 =  x_N;

end
toc
save bin_sec bin_x;
t = int8(0);
key = zeros(n/4,1,'uint8');
for ind1 = 1 : n/4

    for ind2 = 1 : 4
    key(ind1) = key(ind1) + bin_x(ind2*ind1)* 2 ^ (ind2-1);
    end

end


text:
function a = dif2_inverse ( n )

%*****************************************************************************80
%
%% DIF2_INVERSE returns the inverse of the DIF2 matrix.
%
%  Formula:
%
%    if ( I <= J )
%      A(I,J) = I * (N-J+1) / (N+1)
%    else
%      A(I,J) = J * (N-I+1) / (N+1)
%
%  Example:
%
%    N = 4
%
%            4 3 2 1
%    (1/5) * 3 6 4 2
%            2 4 6 3
%            1 2 3 4
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    det ( A ) = 1 / ( N + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      if ( i <= j )
        a(i,j) = ( i * ( n - j + 1 ) ) / ( n + 1 );
      else
        a(i,j) = ( j * ( n - i + 1 ) ) / ( n + 1 );
      end
    end
  end

  return
end


text:
function [cmap] = gat(nColors, cZero, useWhite)
if nargin < 1
    nColors = 255;
end
if nargin < 2
    cZero = 1;
end
if nargin < 3
    useWhite  = 1;
end

clim = caxis();
if (clim(1)>0 || clim(2)<0)
    cZero = 0;
end
if cZero
    center = round(clim(2)*nColors/(clim(2)-clim(1)));
    N = 2*center;
else
    N = nColors;
    center = round(nColors/2);
end

if nColors == 0
    cmap = [];
    return
end

try
    L = ones(nColors,1)*100;
    a = ones(nColors,1)*0;
    b = ones(nColors,1)*0;

    % RED     100  128  128
    % VIOLET  100  128 -128
    % GREEN     x -128  128
    % BLUE      x -128 -128

    % BLACK => RED => GREEN =>  BLUE => VIOLET

    % BLACK to RED
    nVal = floor(N*5/32);
    Lab1 = [20 20 10];
    Lab2 = [40 100 128];
    step = 1/(nVal-1);
    L(1:nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
    a(1:nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
    b(1:nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
    s = nVal;

    % RED to YELLOW
    nVal = floor(N*6/32);
    Lab1 = [40 100 120];
    Lab2 = [100 -20 128];
    step = 1/(nVal);
    L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','pchip');
    a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
    b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
    s = s+nVal;

    % YELLOW to GREEN
    nVal = floor(N*3/32);
    Lab1 = [100 0 128];
    Lab2 = [75 -80 80];
    step = 1/(nVal);
    L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','pchip');
    a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
    b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
    s = s+nVal;

    % GREEN to useWhite
    if (useWhite)
        nVal = floor(N*2/32);
        Lab1 = [75 -80 80];
        Lab2 = [100 0 0];
        step = 1/(center-s);
        L(s:center) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:center) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:center) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
        s = center;
    end

    % HALF BAR ------- ------- ------- ------- ------- ------- ------- -------

    % GREEN to BLUE
    if (~useWhite)
        nVal = floor(N*3/32);
        Lab1 = [75 -80 80];
        Lab2 = [100 -80  0];
        step = 1/(nVal);
        L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
        s = s+nVal;
    end

    if cZero
        N = 2*(nColors-center);
    end

    if N > 0
        if (~useWhite)
            nVal = floor(N*2/32);
            Lab1 = [100 -80 0];
            Lab2 = [80 40 -128];
            step = 1/(nVal);
            L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
            a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
            b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
            s = s+nVal;
        end
        % HALF BAR ------- ------- ------- ------- ------- ------- ------- -------

        if (useWhite)
            % useWhite to BLUE
            nVal = floor(N*2/32);
            Lab1 = [100 0 0];
            Lab2 = [80 40 -128];
            step = 1/(nVal);
            L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
            a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
            b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
            s = s+nVal;
        end

        % BLUE to PINK
        nVal = floor(N*3/32);
        Lab1 = [80 40 -128];
        Lab2 = [50 -80 -80];
        step = 1/(nVal);
        L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
        s = s+nVal;

        % BLUE to VIOLET
        nVal = floor(N*4/32);
        Lab1 = [50 -80 -80];
        Lab2 = [0 100 -128];
        step = 1/(nVal);
        L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
        s = s+nVal;

        % BLUE to VIOLET
        nVal = floor(N*4/32);
        Lab1 = [0  100 -128];
        Lab2 = [40 128 -128];
        step = 1/(nVal);
        L(s:s+nVal) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:s+nVal) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:s+nVal) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
        s = s+nVal;

        % BLUE to VIOLET
        nVal = floor(N*4/32);
        Lab1 = [40 128 -128];
        Lab2 = [25 10 -20];
        step = 1/(nColors-s);
        L(s:end) = interp1([0 1]', [Lab1(1) Lab2(1)], [0:step:1]','linear');
        a(s:end) = interp1([0 1]', [Lab1(2) Lab2(2)], [0:step:1]','linear');
        b(s:end) = interp1([0 1]', [Lab1(3) Lab2(3)], [0:step:1]','linear');
    end

    %L = [[0:100/(nColors/2-1):100]' ; flipud([0:100/(nColors/2-1):100]')];
    %L = flipud([0:100/(nColors-1):100]');
    cmap = flipud(Lab2RGB(L, a, b )./255);
catch
    cmap = jet(1024);

end
end

function [R, G, B] = Lab2RGB(L, a, b)

var_Y = ( L + 16 ) ./ 116;
var_X = a ./ 500 + var_Y;
var_Z = var_Y - b ./ 200;

pos = (var_Y.^3 > 0.008856);
var_Y(pos) = var_Y(pos).^3;
var_Y(~pos) = (var_Y(~pos) - 16 / 116 ) ./ 7.787;

pos = (var_X.^3 > 0.008856);
var_X(pos) = var_X(pos).^3;
var_X(~pos) = (var_X(~pos) - 16 / 116 ) ./ 7.787;

pos = (var_Z.^3 > 0.008856);
var_Z(pos) = var_Z(pos).^3;
var_Z(~pos) = (var_Z(~pos) - 16 / 116 ) ./ 7.787;

% Observer= 2?, Illuminant= D65
ref_X =  95.047;
ref_Y = 100.000;
ref_Z = 108.883;

X = ref_X .* var_X;
Y = ref_Y .* var_Y;
Z = ref_Z .* var_Z;

var_X = X ./ 100;        % X from 0 to  95.047      (Observer = 2?, Illuminant = D65)
var_Y = Y ./ 100;        % Y from 0 to 100.000
var_Z = Z ./ 100;        % Z from 0 to 108.883

var_R = var_X .*  3.2406 + var_Y .* -1.5372 + var_Z .* -0.4986;
var_G = var_X .* -0.9689 + var_Y .*  1.8758 + var_Z .*  0.0415;
var_B = var_X .*  0.0557 + var_Y .* -0.2040 + var_Z .*  1.0570;

pos = ( var_R > 0.0031308 );
var_R(pos) = 1.055 .* ( var_R(pos) .^ ( 1 / 2.4 ) ) - 0.055;
var_R(~pos) = 12.92 .* var_R(~pos);

pos = ( var_G > 0.0031308 );
var_G(pos) = 1.055 .* ( var_G(pos) .^ ( 1 / 2.4 ) ) - 0.055;
var_G(~pos) = 12.92 .* var_G(~pos);

pos = ( var_B > 0.0031308 );
var_B(pos) = 1.055 .* ( var_B(pos) .^ ( 1 / 2.4 ) ) - 0.055;
var_B(~pos) = 12.92 .* var_B(~pos);

R = max(0,min(var_R .* 255,255));
G = max(0,min(var_G .* 255,255));
B = max(0,min(var_B .* 255,255));

if ((nargout == 1) || (nargout == 0))
    R = [R,G,B];
end

end

author: goGPS-Project