text:
function [calculation_times]=PIVlab_commandline_parallel_speed(directory,cpu_core_number)

nr_of_cores =cpu_core_number; % integer, 1 means single core, greater than 1 means parallel
if nr_of_cores > 1
	try
		local_cluster=parcluster('local'); % single node
		corenum =  local_cluster.NumWorkers ; % fix : get the number of cores available
	catch
		warning('on');
		warning('parallel local cluster can not be created, assigning number of cores to 1');
		nr_of_cores = 1;
	end
end

if nr_of_cores > 1
	if pivparpool('size')~=nr_of_cores
		pivparpool('close');
		pivparpool('open',nr_of_cores);
	end
end

%% Create list of images inside user specified directory

% default directory: PIVlab/Examples
addpath(fileparts(mfilename('fullpath')));
%disp(fileparts(mfilename('fullpath')))

suffix='*.jpg'; %*.bmp or *.tif or *.jpg or *.tiff or *.jpeg
%disp(['Looking for ' suffix ' files in the selected directory.']);
direc = dir([directory,filesep,suffix]); filenames={};
[filenames{1:length(direc),1}] = deal(direc.name);
filenames = sortrows(filenames); %sort all image files
amount = length(filenames);

%% Standard PIV Settings
s = cell(15,2); % To make it more readable, let's create a "settings table"
%Parameter                          %Setting           %Options
s{1,1}= 'Int. area 1';              s{1,2}=128;         % window size of first pass
s{2,1}= 'Step size 1';              s{2,2}=64;         % step of first pass
s{3,1}= 'Subpix. finder';           s{3,2}=1;          % 1 = 3point Gauss, 2 = 2D Gauss
s{4,1}= 'Mask';                     s{4,2}=[];         % If needed, generate via: imagesc(image); [temp,Mask{1,1},Mask{1,2}]=roipoly;
s{5,1}= 'ROI';                      s{5,2}=[];         % Region of interest: [x,y,width,height] in pixels, may be left empty
s{6,1}= 'Nr. of passes';            s{6,2}=2;          % 1-4 nr. of passes
s{7,1}= 'Int. area 2';              s{7,2}=64;         % second pass window size
s{8,1}= 'Int. area 3';              s{8,2}=16;         % third pass window size
s{9,1}= 'Int. area 4';              s{9,2}=16;         % fourth pass window size
s{10,1}='Window deformation';       s{10,2}='*linear'; % '*spline' is more accurate, but slower
s{11,1}='Repeated Correlation';     s{11,2}=0;         % 0 or 1 : Repeat the correlation four times and multiply the correlation matrices.
s{12,1}='Disable Autocorrelation';  s{12,2}=0;         % 0 or 1 : Disable Autocorrelation in the first pass.
s{13,1}='Correlation style';        s{13,2}=0;         % 0 or 1 : Use circular correlation (0) or linear correlation (1).
s{14,1}='Repeat last pass';   s{14,2}=0; % 0 or 1 : Repeat the last pass of a multipass analyis
s{15,1}='Last pass quality slope';   s{15,2}=0.025; % Repetitions of last pass will stop when the average difference to the previous pass is less than this number.

%% Standard image preprocessing settings
p = cell(10,1);
%Parameter                       %Setting           %Options
p{1,1}= 'ROI';                   p{1,2}=s{5,2};     % same as in PIV settings
p{2,1}= 'CLAHE';                 p{2,2}=1;          % 1 = enable CLAHE (contrast enhancement), 0 = disable
p{3,1}= 'CLAHE size';            p{3,2}=50;         % CLAHE window size
p{4,1}= 'Highpass';              p{4,2}=0;          % 1 = enable highpass, 0 = disable
p{5,1}= 'Highpass size';         p{5,2}=15;         % highpass size
p{6,1}= 'Clipping';              p{6,2}=0;          % 1 = enable clipping, 0 = disable
p{7,1}= 'Wiener';                p{7,2}=0;          % 1 = enable Wiener2 adaptive denoise filter, 0 = disable
p{8,1}= 'Wiener size';           p{8,2}=3;          % Wiener2 window size
p{9,1}= 'Minimum intensity';     p{9,2}=0.0;        % Minimum intensity of input image (0 = no change)
p{10,1}='Maximum intensity';     p{10,2}=1.0;       % Maximum intensity on input image (1 = no change)

%% other settings
pairwise = 1; % 0 for [A+B], [B+C], [C+D]... sequencing style, and 1 for [A+B], [C+D], [E+F]... sequencing style

%% PIV analysis loop
if pairwise == 1
	if mod(amount,2) == 1 %Uneven number of images?
		disp('Image folder should contain an even number of images.')
		%remove last image from list
		amount=amount-1;
		filenames(size(filenames,1))=[];
	end

	%disp(['Found ' num2str(amount) ' images (' num2str(amount/2) ' image pairs).'])
	x=cell(amount/2,1);
	y=x;
	u=x;
	v=x;
else
	%disp(['Found ' num2str(amount) ' images'])
	x=cell(amount-1,1);
	y=x;
	u=x;
	v=x;
end

typevector=x; %typevector will be 1 for regular vectors, 0 for masked areas
correlation_map=x; % correlation coefficient

%% Pre-load the image names out side of the parallelizable loop
slicedfilename1=cell(0);
slicedfilename2=cell(0);
j = 1;
for i=1:1+pairwise:amount-1
	slicedfilename1{j}=filenames{i}; % begin
	slicedfilename2{j}=filenames{i+1}; % end
	j = j+1;
end

dummy=imread(fullfile(directory,filenames{1}));
disp(['    Current image size:    ' num2str(size(dummy,1)) ' px'])

tstart=tic;
%% Main PIV analysis loop:
% parallel
disp('Please wait....')
if nr_of_cores > 1

	parfor i=1:size(slicedfilename1,2)  % index must increment by 1

		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,false);
	end
else % sequential loop

	for i=1:size(slicedfilename1,2)  % index must increment by 1

		[x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
			piv_analysis(directory, slicedfilename1{i}, slicedfilename2{i},p,s,nr_of_cores,false);

		%disp([int2str((i+1)/amount*100) ' %']);

	end
end

%% PIV postprocessing loop
% Standard image post processing settings

r = cell(6,1);
%Parameter     %Setting                                     %Options
r{1,1}= 'Calibration factor, 1 for uncalibrated data';      r{1,2}=1;                   % Calibration factor for u
r{2,1}= 'Calibration factor, 1 for uncalibrated data';      r{2,2}=1;                   % Calibration factor for v
r{3,1}= 'Valid velocities [u_min; u_max; v_min; v_max]';    r{3,2}=[-50; 50; -50; 50];  % Maximum allowed velocities, for uncalibrated data: maximum displacement in pixels
r{4,1}= 'Stdev check?';                                     r{4,2}=1;                   % 1 = enable global standard deviation test
r{5,1}= 'Stdev threshold';                                  r{5,2}=7;                   % Threshold for the stdev test
r{6,1}= 'Local median check?';                              r{6,2}=1;                   % 1 = enable local median test
r{7,1}= 'Local median threshold';                           r{7,2}=3;                   % Threshold for the local median test

u_filt=cell(size(u));
v_filt=cell(size(v));
typevector_filt=typevector;

if nr_of_cores >1 % parallel

	if pivparpool('size')~=nr_of_cores
		pivparpool('open',nr_of_cores);
	end

	parfor PIVresult=1:size(x,1)

		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);

	end

else % sequential loop

	for PIVresult=1:size(x,1)

		[u_filt{PIVresult,1}, v_filt{PIVresult,1},typevector_filt{PIVresult,1}]= ...
			post_proc_wrapper(u{PIVresult,1},v{PIVresult,1},typevector{PIVresult,1},r,true);

	end

end
calculation_times=toc(tstart);
%% clean up parallel pool, and cluster
%{
if nr_of_cores >1 % parallel
	poolobj = gcp('nocreate'); % GET the current parallel pool
	if ~isempty(poolobj ); delete(poolobj );end
	clear local_cluster;
end
%}
%%
%save(fullfile(directory, [filenames{1} '_' filenames{end} '_' num2str(amount) '_frames_result_.mat']));

%%
%clearvars -except p s r x y u v typevector directory filenames u_filt v_filt typevector_filt correlation_map
disp('DONE.')

end
function [u_filt, v_filt,typevector_filt] = post_proc_wrapper(u,v,typevector,post_proc_setting,paint_nan)
% wrapper function for PIVlab_postproc

% INPUT
% u, v: u and v components of vector fields
% typevector: type vector
% post_proc_setting: post processing setting
% paint_nan: bool, whether to interpolate missing data

% OUTPUT
% u_filt, v_filt: post-processed u and v components of vector fields
% typevector_filt: post-processed type vector

[u_filt,v_filt] = PIVlab_postproc(u,v, ...
	post_proc_setting{1,2},...
	post_proc_setting{2,2},...
	post_proc_setting{3,2},...
	post_proc_setting{4,2},...
	post_proc_setting{5,2},...
	post_proc_setting{6,2},...
	post_proc_setting{7,2});

typevector_filt = typevector; % initiate
typevector_filt(isnan(u_filt))=2;
typevector_filt(isnan(v_filt))=2;
typevector_filt(typevector==0)=0; %restores typevector for mask

if paint_nan
	u_filt=inpaint_nans(u_filt,4);
	v_filt=inpaint_nans(v_filt,4);
end

end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function [f,g,H] = autoHess(x,type,funObj,varargin)
% Numerically compute Hessian of objective function from gradient values

p = length(x);

if type == 1
	% Use finite differencing
	mu = 2*sqrt(1e-12)*(1+norm(x));

	[f,g] = funObj(x,varargin{:});
	diff = zeros(p);
	for j = 1:p
		e_j = zeros(p,1);
		e_j(j) = 1;
		[f diff(:,j)] = funObj(x + mu*e_j,varargin{:});
	end
	H = (diff-repmat(g,[1 p]))/mu;
elseif type == 3 % Use Complex Differentials
	mu = 1e-150;

	diff = zeros(p);
	for j = 1:p
		e_j = zeros(p,1);
		e_j(j) = 1;
		[f(j) diff(:,j)] = funObj(x + mu*i*e_j,varargin{:});
	end
	f = mean(real(f));
	g = mean(real(diff),2);
	H = imag(diff)/mu;
else % Use central differencing
	mu = 2*sqrt(1e-12)*(1+norm(x));

	f1 = zeros(p,1);
	f2 = zeros(p,1);
	diff1 = zeros(p);
	diff2 = zeros(p);
	for j = 1:p
		e_j = zeros(p,1);
		e_j(j) = 1;
		[f1(j) diff1(:,j)] = funObj(x + mu*e_j,varargin{:});
		[f2(j) diff2(:,j)] = funObj(x - mu*e_j,varargin{:});
	end
	f = mean([f1;f2]);
	g = mean([diff1 diff2],2);
	H = (diff1-diff2)/(2*mu);
end

% Make sure H is symmetric
H = (H+H')/2;

if 0 % DEBUG CODE
	[fReal gReal HReal] = funObj(x,varargin{:});
	[fReal f]
	[gReal g]
	[HReal H]
	pause;
end

author: zouchuhang

repo: LayoutNet


text:
function rgb = oklab2rgb(ok)
  % OKLAB2RGB Convert colors stored in OKLab color space to RGB color space
  %
  % rgb = oklab2rgb(ok)
  %
  % Inputs:
  %   ok  #h by #w by 3 image of colors or #ok by 3 list of colors
  % Outputs:
  %   rgb  converted ok colors (same size as ok)
  %
  % See also: rgb2lin, lin2rgb, oklab2lin , rgb2oklab, lin2oklab
  lin = oklab2lin(ok);
  rgb = lin2rgb(lin);
end

author: alecjacobson

repo: gptoolbox


text:
function computeAllPredictionTime_HN(pathModels,pathResults,training,param,maxOrder,nBoot,seed)
% -------------------------------------------------------------------------
% function computeAllPrediction_HN(pathModels,pathResults,training,param,maxOrder,nBoot,imbalance,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for a given
% feature set type, and for all model orders of all experiments with
% different degrees of freedom. See ref. [1,2] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 2. pathModels: Full path to the Models folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/MODELS'
% 2. pathResults: Full path to the Results folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/COHORT-BASED-RESULTS/Experiment1/RESULTS'
% 3. training: Structure defining all parameters for the given experiment
%              to perform. See masterScript_HN.m for more details.
% 4. param: Cell of two strings, defining 1) The feature set type/name; and
%           2) the outcome to model
% 5. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 6. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 7. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 8. batchNum: (optional input). If present, integer that specifies the
%              batch number for parallelization purposes.
%              --> Ex: 6
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named 'RESULTS' in the
% corresponding folder of each experiment (e.g. 'Experiment1', 'Experiment2', etc.)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off

cd(pathModels)
fSetName = param{1};
outcomeName = param{2};
outcome = training.(outcomeName).timeToEvent;
censoring = 1 - training.(outcomeName).outcome;
models = load(['MODELS_',fSetName,'_',outcomeName]); models = struct2cell(models); models = models{1};

tic
fprintf(['\n --> COMPUTING PREDICTION PERFORMANCE (MODEL ORDERS OF 1 to % u) FOR  "',fSetName,'" FEATURE SET AND "',outcomeName,'" OUTCOME ... '],maxOrder)
for j = 1:maxOrder
    orderName = ['Order',num2str(j)];
    data = models.(orderName).Data;
    [orderResults] = predictionPerformanceEstimationTime_HN(data,outcome,nBoot,censoring,seed);
    results.(orderName) = orderResults;
    results.(orderName).Data = models.(orderName).Data;
    results.(orderName).Name = models.(orderName).Name;
end
cd(pathResults), save(['RESULTS_',fSetName,'_',outcomeName],'results')
fprintf('DONE!\n'), toc

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% This function selects events around a seismic
% station that fullfill certain criteria;

report_this_filefun(mfilename('fullpath'));

load /Seis/obelix/stefan/split_data2/ak8897.mat
ty1 = '+'
ty2 = 'o'
ty2 = 'x'
par1 = 100;minmag = 4;

los = input('Longitude: ')
las = input('Latitude: ')

pl = plot(los,las,'rs')
set(pl,'LineWidth',1.0,'MarkerSize',10,...
    'MarkerFaceColor','y','MarkerEdgeColor','k');

l = sqrt(((a.Longitude-los)*cos(pi/180*las)*111).^2 + ((a.Latitude-las)*111).^2) ;
l2 = a.Magnitude >= 0.0 & a.Depth >= l;
%l2 = a.Magnitude >=2.0 & a.Depth <= 30 & l < 100;
a = a.subset(l2);
update(mainmap())
pl = plot(los,las,'rs')
set(pl,'LineWidth',1.0,'MarkerSize',10,...
    'MarkerFaceColor','r','MarkerEdgeColor','y');

author: CelsoReyes


text:
function value = i4_factorial2 ( n )

%*****************************************************************************80
%
%% I4_FACTORIAL2 computes the factorial N!!
%
%  Formula:
%
%    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
%                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the double factorial function.
%    If N is less than 1, VALUE is returned as 1.
%
%    Output, integer VALUE, the value of N!!.
%
  value = 1;

  if ( n < 1 )
    return
  end

  while ( 1 < n )
    value = value * n;
    n = n - 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QPSK demonstration packet-based transceiver for Chilipepper
% Toyon Research Corp.
% http://www.toyon.com/chilipepper.php
% Created 10/17/2012
% embedded@toyon.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of a Costas Loop. Refer to:
% Telecommunications Breakdown: Concepts of Communication Transmitted via
% Software-Defined Radio C. Richard Johnson
% We employ a hard-decision feedback in order to get rid of the loop
% filters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#codegen
function [z_i_out, z_q_out, fe] = qpsk_rx_foc(y_i, y_q, mu_in)

    persistent phi

    lSin = SIN;
    lCos = COS;

    if isempty(phi)
        phi = 0;
    end

    mu = mu_in/2^12;

    % create the VCO signal
    if phi >= 1
        phi = phi - 1;
    end
    if phi < 0
        phi = phi + 1;
    end

    phi12 = round(phi*2^12)+1;
    if phi12 >= 2^12
        phi12 = 1;
    end
    if phi12 < 0
        phi12 = 0;
    end
    f_i = lCos(phi12+1);
    f_q = lSin(phi12+1);
    ti1 = y_i*f_i;
    ti2 = y_q*f_q;
    tq1 = y_q*f_i;
    tq2 = -y_i*f_q;
    z_i = ti1 + ti2;
    z_q = tq1 + tq2;

    % generate the error term to drive VCO generateion
    if z_q < 0
        tf = -z_i;
    else
        tf = z_i;
    end
    if z_i < 0
        bf = -z_q;
    else
        bf = z_q;
    end
    % using sign of error in order to make it gain invariant
    time_diff = tf-bf;
    if time_diff < 0
        e = -1;
    else
        e = 1;
    end

    c = mu*e;
    phiNew = phi - c;
    phi = phiNew;

    fe = phiNew;

    z_i_out = z_i;
    z_q_out = z_q;
end

author: Sable

repo: mcbench-benchmarks


text:
function [delta,h,alpha] = iswnbr(w,thetaSQR)
% [delta,h,alpha] = iswnbr(vSQR,thetaSQR)
%
% ISWNBR  Checks feasibility w.r.t. wide region/neighborhood of Sturm-Zhang.
%   vTAR:= (1-alpha)*max(h,v) projection v onto theta-central region
%   delta = (sqrt(n)/theta) * norm(vTAR - v) / norm(v)
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

% ----------------------------------------
% r = n/thetaSQR
% hSQR = sumwNT/(r-|T|),  hubSQR = sumwNT/(r-|T| - |Q|)
% sumdifv = h*|T| - sumvT   (sumvT = sum(v_T), growing )
% sumdifw = hSQR*|T| - sumwT
% alpha = sumdifv / (r*h)
% deltaSQR = r * ( 2*alpha-alpha^2 - (1-alpha)^2 * sumdifw/gap )
% WE UPDATE sumdifv AND sumdifw IN A STABLE WAY
% ----------------------------------------
n = length(w); gap = sum(w);
sumwNT = gap;
r = n / thetaSQR;
cardT = 0; wQ = []; sumdifv = 0; sumdifw = 0;
cardQ = n;
hSQR = sumwNT / (r - cardT);  hubSQR = sumwNT / (r-(n-1));
for j = 1:n
    wj = w(j);
    if wj >= hubSQR              % wj >= hubSQR ==> not in T
        cardQ = cardQ - 1;
        hubSQR = sumwNT / (r-cardT-cardQ);
    elseif wj < hSQR             % wj < hSQR ==> in T
        cardT = cardT + 1;
        cardQ = cardQ - 1;
        hubSQR = (1-wj/sumwNT) * hubSQR;
        sumwNT = sumwNT - wj;
        oldhSQR = hSQR;
        hSQR = sumwNT / (r - cardT);
        sumdifw = sumdifw + (oldhSQR-wj) + cardT * (hSQR-oldhSQR);
        sumdifv = sumdifv + (sqrt(oldhSQR)-sqrt(wj)) + ...
            cardT * (sqrt(hSQR)-sqrt(oldhSQR));
    else                    % Inconclusive: j in Q
        wQ = [wQ;wj]; %#ok
    end % if
end % for
% ----------------------------------------
% The same treatment for the Q set, but we
% sort the (presumably short) wQ first.
% ----------------------------------------
if ~isempty(wQ)
    sort(wQ);
    STOP = 0; j = 1;
    while ~STOP
        wj = wQ(j);
        if wj >= hSQR
            STOP = 1;
        else
            cardT = cardT + 1;
            sumwNT = sumwNT - wj;
            oldhSQR = hSQR;
            hSQR = sumwNT / (r - cardT);
            sumdifw = sumdifw + (oldhSQR-wj) + cardT * (hSQR-oldhSQR);
            sumdifv = sumdifv + (sqrt(oldhSQR)-sqrt(wj)) + ...
                cardT * (sqrt(hSQR)-sqrt(oldhSQR));
            j = j+1;
            if j > length(wQ)
                STOP = 1;
            end
        end
    end
end % treatment Q
% ----------------------------------------
% alpha = sumdifv/(r*h)
% deltaSQR = r * ( 2*alpha-alpha^2 - (1-alpha)^2 * sumdifw/gap )
%  (THE ABOVE DIFFERENCE SHOULD NOT BE NUMERICALLY DANGEROUS,
%    SINCE alpha IS *SIGNIF* BIGGER THAN sumdifw/gap )
% ----------------------------------------
h = sqrt(hSQR);
alpha = sumdifv/ (r*h);
deltaSQR = alpha*(2-alpha) - (1-alpha)^2 * sumdifw/gap;
delta = sqrt(r*deltaSQR);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function G = repmat(f, M, N)
%REPMAT   Replicate and tile a CHEBFUN.
%   REPMAT(F, M, N) or REPMAT(F, [M, N]) creates an array-valued CHEBFUN by
%   tiling copies of F. If F is a column CHEBFUN, then REPMAT(F, 1, N) returns
%   an array-valued CHEBFUN with N*SIZE(F, 2) CHEBFUN columns. If F is a row
%   CHEBFUN, REPMAT(F, M, 1) returns an array-valued CHEBFUN with M*size(F, 1).
%
% See also HORZCAT, VERTCAT, CAT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Parse inputs:
if ( nargin == 2 )
    if ( length(M) ~= 2 )
        error('CHEBFUN:CHEBFUN:repmat:notEnoughInputs', ...
            'Requires REPMAT(F, M, N) or REPMAT(F, [M, N]).')
    end
    N = M(2);
    M = M(1);
end

if ( f(1).isTransposed )
    % REPMAT a row CHEBFUN:
    if ( N ~= 1 )
        error('CHEBFUN:CHEBFUN:repmat:row',...
            'Use REPMAT(F, M, 1) to replicate and tile row CHEBFUN objects.')
    else
        G = repmat({f}, M, 1);
        G = vertcat(G{:});
    end
else
    % REPMAT a column CHEBFUN:
    if ( M ~= 1 )
        error('CHEBFUN:CHEBFUN:repmat:col',...
            'Use REPMAT(F, 1, N) to replicate and tile column CHEBFUN objects.')
    else
        G = repmat({f}, 1, N);
        G = horzcat(G{:});
    end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [output, Greg] = dftregistration_min_max_3d(buf1ft,buf2ft,usfac,min_shift,max_shift,phase_flag)
% function [output Greg] = dftregistration(buf1ft,buf2ft,usfac);
% Efficient subpixel image registration by crosscorrelation. This code
% gives the same precision as the FFT upsampled cross correlation in a
% small fraction of the computation time and with reduced memory
% requirements. It obtains an initial estimate of the crosscorrelation peak
% by an FFT and then refines the shift estimation by upsampling the DFT
% only in a small neighborhood of that estimate by means of a
% matrix-multiply DFT. With this procedure all the image points are used to
% compute the upsampled crosscorrelation.
% Manuel Guizar - Dec 13, 2007
%
% Rewrote all code not authored by either Manuel Guizar or Jim Fienup
% Manuel Guizar - May 13, 2016
%
% Modified by Eftychios A. Pnevmatikakis to include upper bound on possible
% shifts - November 1, 2016
%
% Citation for this algorithm:
% Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
% "Efficient subpixel image registration algorithms," Opt. Lett. 33,
% 156-158 (2008).
%
% Inputs
% buf1ft    Fourier transform of reference image,
%           DC in (1,1)   [DO NOT FFTSHIFT]
% buf2ft    Fourier transform of image to register,
%           DC in (1,1) [DO NOT FFTSHIFT]
% usfac     Upsampling factor (integer). Images will be registered to
%           within 1/usfac of a pixel. For example usfac = 20 means the
%           images will be registered within 1/20 of a pixel. (default = 1)
% min_shift Minimum shift in each dimension (3x1 vector). (default = -Inf, no min)
% max_shift Maximum shift in each dimension (3x1 vector). (default = Inf, no max)
%
% Outputs
% output =  [error,diffphase,net_row_shift,net_col_shift]
% error     Translation invariant normalized RMS error between f and g
% diffphase     Global phase difference between the two images (should be
%               zero if images are non-negative).
% net_row_shift net_col_shift   Pixel shifts between images
% Greg      (Optional) Fourier transform of registered version of buf2ft,
%           the global phase difference is compensated for.
%
%
% Copyright (c) 2016, Manuel Guizar Sicairos, James R. Fienup, University of Rochester
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the University of Rochester nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

if ~exist('phase_flag','var')
    phase_flag = true;
end

if ~exist('usfac','var')
    usfac = 1;
end

if ~exist('max_shift','var')
    max_shift = Inf(1,3);
end

if ~exist('min_shift','var');
    min_shift = -max_shift;
end

if isscalar(max_shift); max_shift = max_shift*[1,1,1]; end
if isscalar(min_shift); min_shift = min_shift*[1,1,1]; end

[nr,nc,np]=size(buf2ft);
Nr = ifftshift(-fix(nr/2):ceil(nr/2)-1);
Nc = ifftshift(-fix(nc/2):ceil(nc/2)-1);
Np = ifftshift(-fix(np/2):ceil(np/2)-1);

buf_prod = buf1ft.*conj(buf2ft);
if usfac == 0
    % Simple computation of error and phase difference without registration
    CCmax = sum(buf1ft(:).*conj(buf2ft(:)));
    row_shift = 0;
    col_shift = 0;
    pln_shift = 0;
elseif usfac == 1
    % Single pixel registration
    if phase_flag
        buf_prod = buf_prod./abs(buf_prod);
    end
    CC = ifftn(buf_prod);
    CCabs = abs(CC);
    [~,ind] = max(CCabs(:));
    [row_shift, col_shift, pln_shift] = ind2sub([nr,nc,np],ind);
    %[row_shift, col_shift] = find(CCabs == max(CCabs(:)));
    if Nr(row_shift) > max_shift(1) || Nc(col_shift) > max_shift(2) || Np(pln_shift) > max_shift(3) || Nr(row_shift) < min_shift(1) || Nc(col_shift) < min_shift(2) || Np(pln_shift) < min_shift(3)
        CCabs2 = CCabs;
        CCabs2(Nr>max_shift(1),:,:) = 0;
        CCabs2(:,Nc>max_shift(2),:) = 0;
        CCabs2(:,:,Np>max_shift(3)) = 0;
        CCabs2(Nr<min_shift(1),:,:) = 0;
        CCabs2(:,Nc<max_shift(2),:) = 0;
        CCabs2(:,:,Np<max_shift(3)) = 0;
        [~,ind] = max(CCabs2(:));
        [row_shift, col_shift, pln_shift] = ind2sub([nr,nc,np],ind);
        %[row_shift, col_shift] = find(CCabs == max(CCabs2(:)),1,'first');
    end
    CCmax = CC(row_shift,col_shift,pln_shift)*nr*nc*np;
    % Now change shifts so that they represent relative shifts and not indices
    row_shift = Nr(row_shift);
    col_shift = Nc(col_shift);
    pln_shift = Np(pln_shift);
elseif usfac > 1
    % Start with usfac == 2
    buf_pad = FTpad3d(buf_prod,[2*nr,2*nc,2*np]);
    if phase_flag
        buf_pad = buf_pad./(abs(buf_pad)+1e-10);
    end
    CC = ifftn(buf_pad);
    CCabs = abs(CC);
    [~,ind] = max(CCabs(:));
    [row_shift, col_shift, pln_shift] = ind2sub([2*nr,2*nc,2*np],ind);
    %[row_shift, col_shift] = find(CCabs == max(CCabs(:)),1,'first');
    % Now change shifts so that they represent relative shifts and not indices
    Nr2 = ifftshift(-fix(nr):ceil(nr)-1);
    Nc2 = ifftshift(-fix(nc):ceil(nc)-1);
    Np2 = ifftshift(-fix(np):ceil(np)-1);
    if Nr2(row_shift)/2 > max_shift(1) || Nc2(col_shift)/2 > max_shift(2) || Np2(pln_shift)/2 > max_shift(3) || Nr2(row_shift)/2 < min_shift(1) || Nc2(col_shift)/2 < min_shift(2) || Np2(pln_shift)/2 < min_shift(3)
        CCabs2 = CCabs;
        CCabs2(Nr2/2>max_shift(1),:,:) = 0;
        CCabs2(:,Nc2/2>max_shift(2),:) = 0;
        CCabs2(:,:,Np2/2>max_shift(3)) = 0;
        CCabs2(Nr2/2<min_shift(1),:,:) = 0;
        CCabs2(:,Nc2/2<min_shift(2),:) = 0;
        CCabs2(:,:,Np2/2<min_shift(3)) = 0;
        [~,ind] = max(CCabs2(:));
        [row_shift, col_shift, pln_shift] = ind2sub([2*nr,2*nc,2*np],ind);
        %CCabs2 = CCabs(abs(Nr)<=max_shift(1),abs(Nc)<=max_shift(2));
        %[row_shift, col_shift] = find(CCabs == max(CCabs2(:)),1,'first');
    end
    CCmax = CC(row_shift,col_shift,pln_shift)*nr*nc*np;
    row_shift = Nr2(row_shift)/2;
    col_shift = Nc2(col_shift)/2;
    pln_shift = Np2(pln_shift)/2;
    % If upsampling > 2, then refine estimate with matrix multiply DFT
    if usfac > 2,
        %%% DFT computation %%%
        % Initial shift estimate in upsampled grid
        row_shift = round(row_shift*usfac)/usfac;
        col_shift = round(col_shift*usfac)/usfac;
        pln_shift = round(pln_shift*usfac)/usfac;
        dftshift = fix(ceil(usfac*1.5)/2); %% Center of output array at dftshift+1
        % Matrix multiply DFT around the current shift estimate
        CC = conj(dftups3d(buf2ft.*conj(buf1ft),ceil(usfac*1.5),ceil(usfac*1.5),ceil(usfac*1.5),usfac,...
            dftshift-row_shift*usfac,dftshift-col_shift*usfac,dftshift-pln_shift*usfac));
        % Locate maximum and map back to original pixel grid
        CCabs = abs(CC);
        [~,ind] = max(CCabs(:));
        [rloc, cloc, ploc] = ind2sub(size(CC),ind);
        %[rloc, cloc] = find(CCabs == max(CCabs(:)),1,'first');
        CCmax = CC(rloc,cloc,ploc);
        rloc = rloc - dftshift - 1;
        cloc = cloc - dftshift - 1;
        ploc = ploc - dftshift - 1;
        row_shift = row_shift + rloc/usfac;
        col_shift = col_shift + cloc/usfac;
        pln_shift = pln_shift + ploc/usfac;
    end

    % If its only one row or column the shift along that dimension has no
    % effect. Set to zero.
    if nr == 1
        row_shift = 0;
    end
    if nc == 1
        col_shift = 0;
    end
    if np == 1
        pln_shift = 0;
    end
end

if 0
    rg00 = sum(abs(buf1ft(:)).^2);
    rf00 = sum(abs(buf2ft(:)).^2);
    error = 1.0 - abs(CCmax).^2/(rg00*rf00);
    error = sqrt(abs(error));
else
    error = 0;
end
diffphase = angle(CCmax);

output=[error,diffphase,row_shift,col_shift,pln_shift];

% Compute registered version of buf2ft
if (nargout > 1)&&(usfac > 0)
    [Nc,Nr,Np] = meshgrid(Nc,Nr,Np);
    Greg = buf2ft.*exp(1i*2*pi*(-row_shift*Nr/nr-col_shift*Nc/nc-pln_shift*Np/np));
    Greg = Greg*exp(1i*diffphase);
elseif (nargout > 1)&&(usfac == 0)
    Greg = buf2ft*exp(1i*diffphase);
end
return

function out=dftups3d(in,nor,noc,nop,usfac,roff,coff,poff)
% function out=dftups(in,nor,noc,usfac,roff,coff);
% Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just
% a small region.
% usfac         Upsampling factor (default usfac = 1)
% [nor,noc]     Number of pixels in the output upsampled DFT, in
%               units of upsampled pixels (default = size(in))
% roff, coff    Row and column offsets, allow to shift the output array to
%               a region of interest on the DFT (default = 0)
% Recieves DC in upper left corner, image center must be in (1,1)
% Manuel Guizar - Dec 13, 2007
% Modified from dftus, by J.R. Fienup 7/31/06

% This code is intended to provide the same result as if the following
% operations were performed
%   - Embed the array "in" in an array that is usfac times larger in each
%     dimension. ifftshift to bring the center of the image to (1,1).
%   - Take the FFT of the larger array
%   - Extract an [nor, noc] region of the result. Starting with the
%     [roff+1 coff+1] element.

% It achieves this result by computing the DFT in the output array without
% the need to zeropad. Much faster and memory efficient than the
% zero-padded FFT approach if [nor noc] are much smaller than [nr*usfac nc*usfac]

[nr,nc,np]=size(in);
% Set defaults
if exist('roff', 'var')~=1, roff=0;  end
if exist('coff', 'var')~=1, coff=0;  end
if exist('poff', 'var')~=1, poff=0;  end
if exist('usfac','var')~=1, usfac=1; end
if exist('noc',  'var')~=1, noc=nc;  end
if exist('nor',  'var')~=1, nor=nr;  end
if exist('nop',  'var')~=1, nop=np;  end

%tp3d =  @(x,y,z) reshape(kron(y(:),x(:))*z(:)',length(x),length(y),length(z));

% Compute kernels and obtain DFT by matrix products
%kernc=exp((-1i*2*pi/(nc*usfac))*( ifftshift(0:nc-1).' - floor(nc/2) )*( (0:noc-1) - coff ));
kernc=exp((-1i*2*pi/(nc*usfac))*((0:noc-1).' - coff) * (ifftshift(0:nc-1) - floor(nc/2)) );
kernr=exp((-1i*2*pi/(nr*usfac))*((0:nor-1).' - roff) * (ifftshift(0:nr-1) - floor(nr/2)) );
kernp=exp((-1i*2*pi/(np*usfac))*((0:nop-1).' - poff) * (ifftshift(0:np-1) - floor(np/2)) );

%kernc = exp((-1i*2*pi/(nc*usfac))*tp3d(ifftshift(0:nc-1) - floor(nc/2),(0:noc-1)-coff,(0:noc-1)-coff));
%kernr = exp((-1i*2*pi/(nr*usfac))*tp3d((0:nor-1)-roff,ifftshift(0:nr-1) - floor(nr/2),(0:nor-1)-roff));
%kernp = exp((-1i*2*pi/(np*usfac))*tp3d((0:nop-1)-poff,(0:nop-1)-poff,ifftshift(0:np-1) - floor(np/2)));
out = reshape(kernr*reshape(in,nr,[]),nor,nc,np);
out = permute(out,[2,1,3]);
out = reshape(kernc*reshape(out,nc,[]),noc,nor,np);
out = permute(out,[2,1,3]);
out = permute(out,[3,2,1]);
out = reshape(kernp*reshape(out,np,[]),noc,nor,nop);
out = permute(out,[3,2,1]);

%out=kernr*in*kernc;
return

function [ imFTout ] = FTpad3d(imFT,outsize)
% imFTout = FTpad3d(imFT,outsize)
% Pads or crops the Fourier transform to the desired ouput size. Taking
% care that the zero frequency is put in the correct place for the output
% for subsequent FT or IFT. Can be used for Fourier transform based
% interpolation, i.e. dirichlet kernel interpolation.
%
%   Inputs
% imFT      - Input complex array with DC in [1,1,1]
% outsize   - Output size of array [ny nx np]
%
%   Outputs
% imout   - Output complex image with DC in [1,1,1]
% Manuel Guizar - 2014.06.02
% modified by Eftychios A. Pnevmatikakis to include 3d analysis and max
% shifts - 2016.11.03

Nout = outsize;
Nin = size(imFT);
imFT = fftshift(imFT);
center = floor(size(imFT)/2)+1;
imFTout = zeros(outsize);
centerout = floor(size(imFTout)/2)+1;

% imout(centerout(1)+[1:Nin(1)]-center(1),centerout(2)+[1:Nin(2)]-center(2)) ...
%     = imFT;
cenout_cen = centerout - center;

imFTout(max(cenout_cen(1)+1,1):min(cenout_cen(1)+Nin(1),Nout(1)),max(cenout_cen(2)+1,1):min(cenout_cen(2)+Nin(2),Nout(2)),max(cenout_cen(3)+1,1):min(cenout_cen(3)+Nin(3),Nout(3))) ...
    = imFT(max(-cenout_cen(1)+1,1):min(-cenout_cen(1)+Nout(1),Nin(1)),max(-cenout_cen(2)+1,1):min(-cenout_cen(2)+Nout(2),Nin(2)),max(-cenout_cen(3)+1,1):min(-cenout_cen(3)+Nout(3),Nin(3)));

%imFTout2 = padarray(imFT,cenout_cen,0,'pre');
%imFTout = padarray(imFTout2,cenout_cen-mod(size(imFT),2),0,'post');
imFTout = ifftshift(imFTout)*Nout(1)*Nout(2)*Nout(3)/(Nin(1)*Nin(2)*Nin(3));
return

author: flatironinstitute

repo: NoRMCorre


text:
function histogram_discrete_test ( )

%*****************************************************************************80
%
%% HISTOGRAM_DISCRETE_TEST tests the HISTOGRAM_DISCRETE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HISTOGRAM_DISCRETE_TEST:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the HISTOGRAM_DISCRETE library.\n' );

  setup_discrete_test ( );
  pdf_discrete_test ( );
  cdf_discrete_test ( );

  s_num = 50;
  bigger_test ( s_num );

  s_num = 200;
  gaussian_test ( s_num );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HISTOGRAM_DISCRETE_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function main_image = PTKGetMainRegionExcludingBorder(threshold_image, minimum_region_volume_mm3, reporting)
    % PTKGetMainRegionExcludingBorder. Finds the largest connected region in a
    % binary 3D volume, excluding any regions which touch the borders in the
    % first and second dimensions.
    %
    % Syntax:
    %     main_image = PTKGetMainRegionExcludingBorder(threshold_image, reporting)
    %
    % Inputs:
    %     threshold_image - a binary 3D volume as a PTKImage.
    %
    %     minimum_region_volume_mm3 (optional) - ignore any regions below this
    %         threshold
    %
    %     reporting (optional) - an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    % Outputs:
    %     main_image - a PTKImage binary volume of the found region
    %
    %
    %     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.
    %

    if ~isa(threshold_image, 'PTKImage')
        error('Requires a PTKImage as input');
    end

    if nargin < 2 || isempty(minimum_region_volume_mm3)
        minimum_region_volume_mm3 = 0;
    end

    if nargin < 3
        reporting = CoreReportingDefault;
    end

    reporting.ShowProgress('Finding lung region');

    bordered_image = true(threshold_image.ImageSize + [2 2 2]);
    bordered_image(2:end-1, 2:end-1, 2:end-1) = threshold_image.RawImage;
    bordered_image(2:end-1, 2:end-1, 1) = false;
    bordered_image(2:end-1, 2:end-1, end) = false;

    image_opening_params_mm = [0, 2, 4, 6];
    image_opening_index = 1;
    still_searching = true;

    while still_searching

        [results, CC] = OpenAndGetRegions(bordered_image, image_opening_params_mm(image_opening_index), threshold_image.VoxelSize, minimum_region_volume_mm3, reporting);

        if numel(results) > 0 || image_opening_index == numel(image_opening_params_mm)
            still_searching = false;
        else
            image_opening_index = image_opening_index + 1;
        end
    end

    if numel(results) == 0
        main_image = [];
        reporting.Error('PTKGetMainRegionExcludingBorder:NoRegionFound', 'No region could be found');
    end

    if numel(results) > 1
        % If more than one region was found (after excluding boundary-touching
        % regions), then check if they are disconnected left and right lungs
        bounding_boxes = regionprops(CC, 'BoundingBox');
        bb_1 = bounding_boxes(results(1)).BoundingBox;
        bb_1 = bb_1([2,1,3,5,4,6]); % Bounding box is [xyz] but Matlab indices are [yxz]
        bb_2 = bounding_boxes(results(2)).BoundingBox;
        bb_2 = bb_2([2,1,3,5,4,6]); % Bounding box is [xyz] but Matlab indices are [yxz]
        image_centre = round(size(bordered_image)/2);

        image_centre_j = image_centre(2);

        if threshold_image.IsCT
            centre_offset = 0;
        else
            % Allow some leeway for MRI images where the lung boundaries may not be so
            % clear
            centre_offset = 10;
        end

        use_both_regions = false;
        if (bb_1(2) >= image_centre_j - centre_offset) && (bb_2(5) < image_centre_j + centre_offset)
            reporting.LogVerbose('I appear to have found 2 disconnected lungs. I am connecting them.');
            use_both_regions = true;
        end

        if (bb_2(2) >= image_centre_j - centre_offset) && (bb_1(5) < image_centre_j + centre_offset)
            reporting.LogVerbose('I appear to have found 2 disconnected lungs. I am connecting them.');
            use_both_regions = true;
        end
    else
        use_both_regions = false;
    end

    bordered_image(:) = false;
    bordered_image(CC.PixelIdxList{results(1)}) = true;
    if use_both_regions
        bordered_image(CC.PixelIdxList{results(2)}) = true;
    end

    main_image = threshold_image.BlankCopy;
    main_image.ChangeRawImage(bordered_image(2:end-1, 2:end-1, 2:end-1));
end

function [results, CC] = OpenAndGetRegions(bordered_image_input, opening_mm, voxel_size, minimum_region_volume_mm3, reporting)

    bordered_image = bordered_image_input;
    if opening_mm > 0
        ball_element = CoreImageUtilities.CreateBallStructuralElement(voxel_size, opening_mm);
        bordered_image = imopen(bordered_image_input > 0, ball_element);
    end

    % Borders are marked as segmented - this ensures all components
    % touching the border are connected and allows us to eliminate them
    % when extracting the lung

    % Obtain connected component matrix
    CC = bwconncomp(bordered_image, 26);

    % Find largest region
    num_pixels = cellfun(@numel, CC.PixelIdxList);
    [sorted_largest_areas, sorted_largest_areas_indices] = sort(num_pixels, 'descend');

    % Remove regions that are below the volume threshold
    pixel_volume = prod(voxel_size);
    sorted_largest_areas_indices = sorted_largest_areas_indices(pixel_volume*sorted_largest_areas >= minimum_region_volume_mm3);

    result_index = 1;
    results = [];
    index_in_sorted_array = 1;
    while (result_index < 3 && index_in_sorted_array <= length(sorted_largest_areas_indices))
        current_region_being_checked = sorted_largest_areas_indices(index_in_sorted_array);
        bordered_image(:) = false;
        bordered_image(CC.PixelIdxList{current_region_being_checked}) = true;
        if (bordered_image(1) || bordered_image(end))
            % This region is connected to the edge
%             reporting.ShowMessage('PTKGetMainRegionExcludingBorder:LargestROIConnectedToExterior', 'The largest region connected with the edge of the volume. I''m assuming this region is outside the body so choosing the next largest region');
        else
            results(result_index) = current_region_being_checked;
            result_index = result_index + 1;
        end
        index_in_sorted_array = index_in_sorted_array + 1;
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
 function y = mod0(x,b)
%function y = mod0(x,b)
%	the usual mod function returns values in the range [0,b)
%	whereas this mod function returns values in the range [-b/2,b/2)
%	which frequently arises in signal processing problems!
%	Jeff Fessler

y = mod(x + b/2, b) - b/2;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function plot_nsat(solution)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nsol=size(solution,1);
if nsol==0
    error('Solution is empty!!!\n');
end

nsat=zeros(nsol,1);
time=zeros(nsol,1);
j=0;
for i=1:nsol
    if dot(solution(i).pos,solution(i).pos)<=0
        continue;
    end
    [~,sow]=time2gpst(solution(i).time);
    time(j+1)=sow;

    nsat(j+1,:)=solution(i).ns;
    j=j+1;
end

if j==0
    error('Solution is empty!!!\n');
end

if j<nsol
    time(j+1:end,:)=[];
    nsat(j+1:end,:)=[];
end

%% plot
maxsat=max(nsat);
H=get(0,'ScreenSize'); w=600; h=450; x=H(3)/2-w/2; y=H(4)/2-h/2;
figure;set(gcf,'Position',[x y w h]);
plot(time,nsat,':','linewidth',1,'color',[0.5,0.5,0.5]);hold on
plot(time,nsat,'.b','linewidth',2,'Markersize',10);
grid on ;set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',0.5);
xlabel('GPS Time (s)'),ylabel('# of satellite');
axis([time(1) time(end) 0 maxsat+5 ]);

return

repo: GINav


text:
clear all;
close all;
more off;
addpath('tools');

epsilon = 1e-5;

x1 = [1.1 0.9 1]';
x2 = [2.2 1.9]';
z  = [1.3 -0.4]';

% get the analytic Jacobian
[e, A, B] = linearize_pose_landmark_constraint(x1, x2, z);

% check the error vector
e_true = [0.135804 0.014684]';
if (norm(e - e_true) > epsilon)
  disp('Your error function seems to return a wrong value');
  disp('Result of your function'); disp(e);
  disp('True value'); disp(e_true);
else
  disp('The computation of the error vector appears to be correct');
end

% compute it numerically
delta = 1e-6;
scalar = 1 / (2*delta);

% test for x1
ANumeric = zeros(2,3);
for d = 1:3
  curX = x1;
  curX(d) += delta;
  err = linearize_pose_landmark_constraint(curX, x2, z);
  curX = x1;
  curX(d) -= delta;
  err -= linearize_pose_landmark_constraint(curX, x2, z);

  ANumeric(:, d) = scalar * err;
end

diff = ANumeric - A;
if max(max(abs(diff))) > epsilon
  disp('Error in the Jacobian for x1');
  disp('Your analytic Jacobian'); disp(A);
  disp('Numerically computed Jacobian'); disp(ANumeric);
  disp('Difference'); disp(diff);
else
  disp('Jacobian for x1 appears to be correct');
end

% test for x2
BNumeric = zeros(2,2);
for d = 1:2
  curX = x2;
  curX(d) += delta;
  err = linearize_pose_landmark_constraint(x1, curX, z);
  curX = x2;
  curX(d) -= delta;
  err -= linearize_pose_landmark_constraint(x1, curX, z);

  BNumeric(:, d) = scalar * err;
end

diff = BNumeric - B;
if max(max(abs(diff))) > epsilon
  disp('Error in the Jacobian for x2');
  disp('Your analytic Jacobian'); disp(B);
  disp('Numerically computed Jacobian'); disp(BNumeric);
  disp('Difference'); disp(diff);
else
  disp('Jacobian for x2 appears to be correct');
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function vectorSez = rotVectToSEZCoords(rVectorECEF, vectorECEF)
    numVectors = size(rVectorECEF,2);

    kHat = repmat([0;0;1], 1, numVectors);
    zHat = vect_normVector(rVectorECEF);
    eHat = vect_normVector(cross(kHat, zHat));
    sHat = vect_normVector(cross(eHat, zHat));

    rotMat = [reshape(sHat, [1, 3, numVectors]);
              reshape(eHat, [1, 3, numVectors]);
              reshape(zHat, [1, 3, numVectors])];

    vectorECEF = reshape(vectorECEF, [3, 1, numVectors]);
    vectorSez = pagemtimes(rotMat, vectorECEF);
    vectorSez = squeeze(vectorSez);
end
%     kHat = [0;0;1];
%     zHat = normVector(rVectorECEF);
%     eHat = normVector(cross(kHat, zHat));
%     sHat = normVector(cross(eHat, zHat));
%
%     rotMat = [sHat'; eHat'; zHat'];
%
%     vectorSez = rotMat * vectorECEF;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%ROT2 SO(2) Rotation matrix
%
% R = ROT2(THETA) is an SO(2) rotation matrix (2x2) representing a rotation of THETA
% radians.
%
% R = ROT2(THETA, 'deg') as above but THETA is in degrees.
%
% See also SE2, TROT2, ISROT2, TRPLOT2, ROTX, ROTY, ROTZ, SO2.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function R = rot2(t, deg)

    if nargin > 1 && strcmp(deg, 'deg')
        t = t *pi/180;
    end

    ct = cos(t);
    st = sin(t);
    R = [
        ct  -st
        st   ct
        ];

repo: paper-simulation


text:
function au_l_d
%
%  Deletes the global data that has been generated by 'au_l_i'.
%
%  Calling sequence:
%
%    au_l_d
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

global LP_L LP_U

clear global LP_L LP_U


text:
function res = cellfun2(funer, structy)
%An alias for cellfun with UniformOutput set to false
if nargin<2
  error('Not enough arguments to cellfun2\n');
end
res = cellfun(funer, structy, 'UniformOutput', false);

author: quantombone

repo: exemplarsvm


text:
function [ n_data, a, b, fx ] = binomial_values ( n_data )

%*****************************************************************************80
%
%% BINOMIAL_VALUES returns some values of the binomial coefficients.
%
%  Discussion:
%
%    The formula for the binomial coefficient is
%
%      C(N,K) = N! / ( K! * (N-K)! )
%
%    In Mathematica, the function can be evaluated by:
%
%      Binomial[n,k]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer A, B, the arguments of the function.
%
%    Output, integer FX, the value of the function.
%
  n_max = 20;

  a_vec = [ ...
     1,  6,  6,  6, 15, ...
    15, 15, 15, 15, 15, ...
    15, 25, 25, 25, 25, ...
    25, 25, 25, 25, 25  ];

  b_vec = [ ...
     0,  1,  3,  5,  1, ...
     3,  5,  7,  9, 11, ...
    13,  1,  3,  5,  7, ...
     9, 11, 13, 15, 17 ];

  fx_vec = [ ...
           1, ...
           6, ...
          20, ...
           6, ...
          15, ...
         455, ...
        3003, ...
        6435, ...
        5005, ...
        1365, ...
         105, ...
          25, ...
        2300, ...
       53130, ...
      480700, ...
     2042975, ...
     4457400, ...
     5200300, ...
     3268760, ...
     1081575 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0;
    b = 0;
    fx = 0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function h = ne(f, g)
%~=   Not equal operator for CHEBFUN objects.
%   H = F ~= G, where F and/or G are CHEBFUN objects, constructs a logical
%   CHEBFUN H which is true (i.e., takes the value 1) where F ~= G, and false
%   (0) elsewhere.
%
% See also EQ, ROOTS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Trivial empty case:
if ( isempty(f) )
    h = f;
    return
elseif ( isempty(g) )
    h = g;
    return
end

% Array-valued?
if ( min(size(f)) > 1 || min(size(g)) > 1 )
	error('CHEBFUN:CHEBFUN:ne:array', ...
        '~= does not support array-valued CHEBFUN objects.');
end

% Call SIGN() to do the work:
h = sign(f - g);

% Get value in interior of each FUN by taking left-sided limit at the breaks:
vals = get(h, 'rval-local');

% Set FUNs that are nonzero to FUNs which are identically 1:
for k = 1:numel(h.funs)
    if ( vals(k) ~= 0 )
        h.funs{k} = 1 + 0*h.funs{k};
    end
end

% pointValues:
ind = isnan(h.pointValues);
h.pointValues = double(logical(h.pointValues(~ind)));

% Tidy the result:
h = merge(h);

end


text:
function dt6FileName = dtiRawFit_Charmed_PD_BPF_highB(dwRaw, bvecs, bvals, outBaseName, bs, fitMethod, adcUnits, xformToAcPc,ROICLmask,BPF,Gmcc,Wmcc,Cmcc,dt6_path)
%
%dt7 = dtiRawFitCharmedwithPDBPF(S_Raw, S_Bvecs, S_Bvals, name, [], 'charmed',[],[],ROICLmask,BPF,G,W,C)%charmed
% dt6FileName = dtiRawFitTensor([dwRaw=uigetfile],
% [bvecsFile=uigetfile], [bvalsFile=uigetfile], [outBaseDir=uigetdir],
% [bootstrapParams=[]], [fitMethod='ls'],
% [adcUnits=dtiGuessAdcUnits], [xformToAcPc=dwRaw.qto_xyz])
%
% Fits a tensor to the raw DW data. The tensors are returned (in [Dxx Dyy
% Dzz Dxy Dxz Dyz] format) and are saved in a dt6 file outBaseName 'dt6.mat'.
%
% If adcUnits is not provided, we'll try to guess based on the magnitude
% of the mean diffusivity. This guess is based on typical values for
% in-vivo human brain tissue. Our preferred units are 'micron^2/msec',
% because they produce very human-friendly numbers (eg. free diffusion of
% water is 3 micron^2/msec). Your adcUnits are determined by the units in
% which your bvals are specified. (adcUnits = 1/bvalUnits) For our GE
% Bammer/Hedehus DTI sequence, the native bvalUnits are sec/mm^2. However,
% we usually divide the bval by 1000, producing msec/micrometer^2 units.
%
% You can also specify a bootstrap to estimate tensor fit uncertainty.
% To skip the bootstrap, make bootstrapParams empty (or specify n=1
% bootstrap samples). bootstrapParams is a struct  with the following
% fields:
%
%   n: number of bootstrap samples. [] or 1 will trigger no bootstrap.
%      We've found that 200-300 gives a good variance estimate, but 500 is
%      probably better.
%
%   nVolsPerRepeat: each bootstrap sample is a permuted dataset. In order
%      to permute the original data in a reasonable way (ie. preserving the
%      # of unique directions and the # of non-DWIs), we need to know the
%      repetion pattern in the input data. For now, we assume that the data
%      are arranged such that the measurements for each repeat are
%      contiguous and that the order of the direction measurements is the
%      same across all the repeats. In that case, one number can specify
%      the approriate pattern for generating the bootstrap permutations.
%      nVolsPerRepeat specifies the number of image volumes to expect per
%      repeat. E.g., if you make 13 measurements per repeat (12 DWIs + 1
%      non DWI), nVolsPerRepeat = 13. Note that the code below will try to
%      deal gracefully with incomplete data for the last repeat.
%
% ===OR===
%
%   permuteMatrix: a cell-array with one entry for each measurement (ie.
%   N = length(permuteMatrix) = size(dwRaw,4). Each entry of this cell
%   array is a 1-d array of indices into the N dwRaw volumes that are
%   valid bootstrap samples of the corresponding volume. E.g:
%      bv = [dlmread(bvecsFile).*repmat(dlmread(bvalsFile),[3 1])];
%      bs.permuteMatrix = {};
%      for(ii=1:size(bv,2))
%         dist1 = sqrt((bv(1,:)-bv(1,ii)).^2+(bv(2,:)-bv(2,ii)).^2+(bv(3,:)-bv(3,ii)).^2);
%         dist2 = sqrt((bv(1,:)+bv(1,ii)).^2+(bv(2,:)+bv(2,ii)).^2+(bv(3,:)+bv(3,ii)).^2);
%         bs.permuteMatrix{ii} = unique([find(dist1<1e-3) find(dist2<1e-3)]);
%      end
%
%   maxMem: The bootstrap tensor fits can go much faster is we do them in
%      large chunks, but not so large that we cause Matlab to use swap
%      space or run out of memory altogether. So set maxMem (specified in
%      either bytes or Gigabytes) to about 80% of the physical RAM that you
%      have available. (defaults to 1GB)
%
%   showProgress: if true, a progress bar will be shown. (defaults to true)
%
% Currently, specifying the bootstrap cases the resulting dt6 file to have
% the following additional variables:
%
%   faStd, mdStd: standard deviations on fa and mean diffusivity.
%
%   pddDisp: dispersion of PDD axes (based on the Watson) in degrees of
%   angle (54 deg is maximum dispersion).
%
%
% fitMethod: the tensor fitting method.
%   'ls': least-squares (default)
%   'me': maximum-entropy method (Dima Kuzmin and Manfred Warmuth, UCSC)
%   'rt': RESTORE robust tensor fitting and outlier rejection:
%         Chang, Jones & Pierpaoli (2005). RESTORE: Robust Estimation of
%         Tensors by Outlier Rejection. Magnetic Resonance in Medicine, v53.
%
% Note that the RESTORE implementation is experimental and needs more
% testing. Also, don't do a bootstrap with RESTORE- that doesn't work yet.
%
% E.g.:
% f = 'raw/dti_g13_b800_aligned.'; out = 'dti06rt';
% dtiRawFitTensor([f 'nii.gz'], [f 'bvecs'], [f 'bvals'], out, [], 'rt');
% % Show outlier count as an overlay on the b0:
% aNi = niftiRead(fullfile(out,'bin','b0.nii.gz'));
% oNi = niftiRead(fullfile(out,'bin','outliers.nii.gz'));
% aIm = mrAnatHistogramClip(double(aNi.data),0.4,0.98);
% oIm = double(sum(oNi.data,4));
% mrAnatOverlayMontage(oIm, oNi.qto_xyz, aIm, aNi.qto_xyz, autumn(256), [1 10], [-34:2:62],[],3,1,false);
%
% TODO:
%
% * add more statistics to bootstrap. We should fit an assymetric
% distribution to the PDD pdf, like the Bingham. Also, we should do
% ksdensity on fa and md and save out a more complete description of the
% PDFs, as they are not well-fit by the normal assumption implicit in the
% standard deviation.
%
% HISTORY:
%
% 2007.03.20 RFD: wrote it.
% 2007.05.30 RFD: added bootstrap option
% 2007.06.02 RFD: cleaned and documented bootstrap code. Seems to work well
% now.
% 2007.06.08 RFD: now save in the new all-NIFTI format. THe dt6 file is now
% just a 'project' file with some notes and the filenames of the actual
% data files.
% 2007.06.14 RFD: NIFTI tensor files weren't respecting the NIFTI-1 spec.
% This is now fixed and the tensor elements are stored in lower-triangular,
% row-wise order (Dxx Dxy Dyy Dxz Dyz Dzz).
% 2007.07.20 AJS: Relative fileanames to the parent directory.
% 2008.09.03 DY & RFD: Implemented 'rt' fitMethod (RESTORE)
% 2008.12.16 DY: Forced useParfor (parallel processing flag) to be false,
% as it seems not to be functional (on Bob's recommendation)
% 2008.12.18 DY & AL: Create summary image of outliers.nii.gz that can be viewed as an image
% when loaded into DTIfiberUI.
% if(license('checkout','distrib_computing_toolbox'))
%     useParfor = true;
% else
%useParfor = false;
% end

if(~exist('dwRaw','var')||isempty(dwRaw))
    [f,p] = uigetfile({'*.nii.gz;*.nii';'*.*'}, 'Select the raw DW NIFTI dataset...');
    if(isnumeric(f)), error('User cancelled.'); end
    dwRaw = fullfile(p,f);
end
if(ischar(dwRaw))
    % dwRaw can be a path to the file or the file itself
    [dataDir,inBaseName] = fileparts(dwRaw);
else
    [dataDir,inBaseName] = fileparts(dwRaw.fname);
end
[junk,inBaseName] = fileparts(inBaseName);
if(isempty(dataDir)), dataDir = pwd; end

if(~exist('bvecs','var')||isempty(bvecs))
    bvecs = fullfile(dataDir,[inBaseName '.bvecs']);
    [f,p] = uigetfile({'*.bvecs';'*.*'},'Select the bvecs file...',bvecs);
    if(isnumeric(f)), disp('User canceled.'); return; end
    bvecs = fullfile(p,f);
end
if(~exist('bvals','var')||isempty(bvals))
    bvals = fullfile(dataDir,[inBaseName '.bvals']);
    [f,p] = uigetfile({'*.bvals';'*.*'},'Select the bvals file...',bvals);
    if(isnumeric(f)), disp('User canceled.'); return; end
    bvals = fullfile(p,f);
end
if(~exist('outBaseName','var')||isempty(outBaseName))
    if(nargout==0)
        outBaseName = fullfile(dataDir,inBaseName);
    else
        outBaseName = [];
    end
end
if(~exist('adcUnits','var'))
    adcUnits = '';
end
if(isempty(outBaseName))
    outBaseName = uigetdir(inBaseName,'Select a directory for the data...');
    if(isnumeric(outBaseName)), disp('User canceled.'); return; end
end
dt6FileName = fullfile(outBaseName, 'dt6.mat');
binDirName = fullfile(outBaseName, 'bin');
if(exist(outBaseName,'dir'))
    if(exist(dt6FileName,'file')||(exist(binDirName,'dir')&&~isempty(dir(fullfile(binDirName,'*.nii*')))))
        q = ['Output dir ' outBaseName ' exists and appears to contain data. Are you sure that you want to overwrite the data files in there?'];
        resp = questdlg(q,'Confirm Overwrite','Yes','Cancel','Cancel');
        if(strcmp(resp,'Cancel')), disp('canceled.'); return; end
        %error(['Output dir ' outBaseName ' exists and appears to contain data- please move it out of the way.']);
        %outBaseName = uigetdir('Select directory for output...',dt6FileName);
        %if(isnumeric(f)), disp('User canceled.'); return; end
        %dt6FileName = fullfile(p,f);
        %[p,f] = fileparts(dt6FileName);
        %binDirName = fullfile(p,f);
    end
end
disp(['data will be saved to ' outBaseName '.']);

if(ischar(dwRaw))
    % dwRaw can be a path to the file or the file itself
    disp(['Loading raw data ' dwRaw '...']);
    dwRaw = niftiRead(dwRaw);
    weLoadedRaw = true;
else
    weLoadedRaw = false;
end

nvols = size(dwRaw.data,4);
if(~exist('xformToAcPc','var')||isempty(xformToAcPc))
    xformToAcPc = dwRaw.qto_xyz;
end

if(~exist('bs','var')||isempty(bs))
    bs.n = 1;
    bs.nVolsPerRepeat = 1;
    % 1 GByte = 2^30
    bs.maxMem = 4*2^30;
    bs.showProgress = false;
else
    if(~isfield(bs,'maxMem')), bs.maxMem = 1*2^30;
    elseif(bs.maxMem<100), bs.maxMem = bs.maxMem*2^30; end
    if(~isfield(bs,'showProgress')), bs.showProgress = true; end
end

if(~exist('fitMethod','var')||isempty(fitMethod))
    fitMethod = 'ls';
end

%% Load the bvecs & bvals
% NOTE: these are assumed to be specified in image space.
% If bvecs are in scanner space, use dtiReorientBvecs and
% dtiRawReorientBvecs.
if(~isnumeric(bvecs))
    %bvecs = dlmread(bvecs, ' ');
    bvecs = dlmread(bvecs);
end
if(~isnumeric(bvals))
    %bvals = dlmread(bvals, ' ');
    bvals = dlmread(bvals);
end

if(size(bvecs,2)~=nvols || size(bvals,2)~=nvols)
    error(['bvecs/bvals: need one entry for each of the ' num2str(nvols) ' volumes.']);
end

%% Get a brain mask
%
disp('Computing brain mask from average b0...');
dwInds = bvals>0;
b0Ims = double(dwRaw.data(:,:,:,~dwInds));
nz = b0Ims>0;
b0Ims(nz) = log(b0Ims(nz));
b0 = exp(mean(b0Ims,4));
clear b0Ims nz;
b0clip = mrAnatHistogramClip(b0,0.4,0.99);
b0 = int16(round(b0));
% We use a liberal brain mask for deciding which tensors to compute, but a
% more conservative mask will be saved so that that the junk outside the
% brain won't be displayed when we view the data.
liberalBrainMask = dtiCleanImageMask(b0clip>0.1&all(dwRaw.data>0,4),10,1,0.25);
liberalBrainMask(all(dwRaw.data==0,4)) = 0;
brainMask = uint8(dtiCleanImageMask(b0clip>0.25));
% make sure the display-purposes brain mask is a subset of the
% tensor-fitting (liberal) brain mask.

brainMask(~liberalBrainMask) = 0;
if(exist('ROICLmask','var') && ~isempty(ROICLmask))
brainMask(ROICLmask)=3;
brainMask(brainMask~=3)=0;
end;
% brainMask(~ROICLmask) = 0;
% brainMask(find(brainMask))=0;
% brainMask(ROICLmask)=1;
clear b0clip badEdgeVox;

%% Reorganize the data to make computations easier
numVols = size(dwRaw.data,4);
brainInds = find(liberalBrainMask);
mask=liberalBrainMask;
if(exist('ROICLmask','var') && ~isempty(ROICLmask))
brainInds = find(brainMask);
mask=brainMask;
end;

data = zeros(numVols,1,length(brainInds));
for(ii=1:numVols)
    tmp = double(dwRaw.data(:,:,:,ii));
    data(ii,1,:) = tmp(brainInds);
end

%% Compute signal noise estimate
%
% According to Henkelman (), the expected signal variance (sigma) can be computed as
% 1.5267 * SD of the background (thermal) noise.
sz = size(dwRaw.data);
x = 10;
y = 10;
z = round(sz(3)/2);
[x,y,z,s] = ndgrid(x-5:x+5,y-5:y:5,z-5:z+5,1:sz(4));
noiseInds = sub2ind(sz,x(:),y(:),z(:),s(:));
sigma = 1.5267 * std(double(dwRaw.data(noiseInds)));

% Memory usage is tight- if we loaded the raw data, clear it now since
% we've made the reorganized copy that we'll use for all subsequent ops.
if(weLoadedRaw), clear dwRaw; end
clear liberalBrainMask;
% Voxels with intensity exactly==0 sometimes occur by chance (usually in a
% artifact region, like a blood vessel) or as a result of eddy/motion
% correction and interpolation. They are a problem because fitting the
% tensor involves taking the log(intensity).
minVal = min(data(data(:)>0));
data(data==0) = minVal;

%% Fit the tensor maps.
%brainMask(~ROICLmask) = 0;

nvox = size(data,3);
% Start with q, which will have a row for each volume, each row having
% three elements: [bvx bvy bvz].
%
% Each row of X corresponds to a DW direction for that volume of the form:
% [1, -b bvx bvx, -b bvy bvy, -b bvz bvz, -2b bvx bvy, -2b bvx bvz, -2b bvy bvz].
%
% The last six values in each row of X: [-bx^2 -by^2 -bz^2 -2bxy -2bxz
% -2byz] are equivalent to equation 1 on p.457 of Basser et al, 2002, NMR
% in Biomedicine. They are the six unique values in the symmetric b-matrix.
%
% Our goal is to use the raw data from each DWI (from dwRaw.data, stored in
% the matlab workspace currently as data) and its corresponding b-matrix
% (currently computed to be X in the matlab workspace) to estimate D using
% multivariate linear regression.

q = (bvecs.*sqrt(repmat(bvals,3,1)))';
X = [ones(numVols,1) -q(:,1).^2 -q(:,2).^2 -q(:,3).^2 -2.*q(:,1).*q(:,2) -2.*q(:,1).*q(:,3) -2.*q(:,2).*q(:,3)];
gof = [];
outliers = [];
switch fitMethod

    case 'charmed'
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%% GET THE Gradient and change there cordinate to speric +innitiate the relevabt Paramer %%%
        g = 42576.0; % kHz/T = (cycles/millisecond)/T
        Par.delta = 41.124;%17.2; % msec (23.1 w/ 40mT = b=1.0), 33.6 for bvals up to 3
        maxG = 50;
        % We need to get the gradient amplitudes that we actually used when
        % scanning. To compute these, we used the default (assumed) value
        % of delta:
        [grads, maxB] = dtiGradsBuildCharmed([], [], [], 0);
        % [grads, maxB, Par.delta, Par.Delta] = dtiGradsBuildCharmed([], [], Par.delta, 0);
        % But we also need the properly-computed Delta
        Par.Delta = maxB/((2*pi*g).^2 * (maxG*1e-9).^2 * Par.delta.^2) + Par.delta/3;
        G = grads*maxG;
        % T/micrometer * cycles/msec/T * msec = cycles/micrometer
        q =  G*1e-9 * g * Par.delta;

        [Q(:,1),Q(:,2) Q(:,3)] =cart2sph(q(1,:)',q(2,:)',q(3,:)');

       % Q(:,1:2)=pi/2-Q(:,1:2);%change tocharmed article convetion

        % Compute the actual bvals
        norm_q = sqrt(sum(q.^2,1));
        bvals = (norm_q.*2*pi).^2*(Par.Delta-Par.delta/3);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% for intiation of some of the fitted parameterswe need to fit the low B value diffution data %%%

BV=unique(sort(bvals));

        normal_b=find(bvals<BV(2)*1.1); % find the law bvals this is the b0 and the secound bvalue as some time there is small digit difference i take all the bval above 10% of the secound bval.
numBV_Vols=length(normal_b);
     q1 = (bvecs(:,normal_b).*sqrt(repmat(bvals(normal_b),3,1)))';
        X1 = [ones(numBV_Vols,1) -q1(:,1).^2 -q1(:,2).^2 -q1(:,3).^2 -2.*q1(:,1).*q1(:,2) -2.*q1(:,1).*q1(:,3) -2.*q1(:,2).*q1(:,3)];

                data1=data(normal_b,:,:);
         Xinv = pinv(X1);
        % The dt values are simply log(data)*Xinv
        % Avoid log of zero warning:
        nz = data1>0;
        data1(nz) = log(data1(nz));
        clear nz;
        A = ndfun('mult', Xinv, data1);

        clear logData;
        %calculate the Tensor for the low b values
        [eigVec,eigVal] = dtiEig(squeeze(A(2:7,1,:))');
        l1 = eigVal(:,1);
        l2 = eigVal(:,2);
        l3 = eigVal(:,3);
        %find the orintation of the firist eigVal and change its cordinate to
        %speric
        eig_L1(:,:)=eigVec(:,[1 2 3],1);
        [Q0(:,1),Q0(:,2) Q0(:,3)] =cart2sph(eig_L1(:,1),eig_L1(:,2),eig_L1(:,3));%
        %Q0(:,1:2)=pi/2-Q0(:,1:2); %change tocharmed article convetion

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

        %%%%% Initiate the fitted parameters

        %the ftted paramters:
        %         theta_H=x(1);  %the  orientation of the hindererd component  from the qradients pitch from z axis
        %         phi_H=x(2);      %the  orientation of the hindererd
        %                           component  from the qradients
        %                       rotation in x-y plane see Assaf et al MRM 2004 p 977-8
        %         f_h=x(3);     the hidered component
        %         Lh_par=x(4);  the parallel Diffusion coefficient in the hindererd component
        %         Lh_per=x(5);  the ratio of the perpendicular ADC to the parallel ADC
        %                       in the hindered component to that of the x(5)=Lh_par/Lh_per;
        %         Dr_par=x(6);  the parallel Diffusion coefficient in the restricted component
        %         N=x(7);       %the noise floor
        %         theta_R=x(8);    %the cylinder orientation in it the diff is
        %                        restricted:pitch from z axis
        %         phi_R=x(9);     %the cylinder orientation in it the diff is restricted:
        %                       rotation in x-y plane see Assaf et al MRM 2004 p 977-8

        options = optimset('LargeScale','on','LevenbergMarquardt','on', 'Display', 'off', 'MaxIter', 100);
        %options = optimset('Display', 'off', 'MaxIter', 100);

        lb=[-10*pi -10*pi   1e-6  1e-6  1e-6  1e-6  1e-6 -10*pi -10*pi ]; %low bondery
        ub=[10*pi   10*pi    1     3     1     3    0.5   10*pi  10*pi ]; %up bondery

        dwInds = bvals>0;
        B0 = mean(data(~dwInds,:),1);
        Dw = data(dwInds,:);
        Par.Q=Q(dwInds,:);
        Par.Dr_per=1; % the restricted perpedicular diffusion coefficient this is a case
        Par.R=0.5;    % the radius of the restricted cilinder. this is a case
        Par.tau=119.4/2; %TE/2 from the charmed data heder file
%         if check==1;
%             %try it on ii=35576 this a CC voxcel
%         normal_b1=find(bvals<BV(2)*1.1 & bvals>0);
%         Dw = data(normal_b1,:);
%         Par.Q=Q(normal_b1,:);
%         end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%fit by charmed model%%%%%
        for(ii=1:nvox)
            raw=Dw(:,ii)/B0(ii);

            x0(1:2)=Q0(ii,1:2); %the orientation of the low b val eigVec1
            x0(8:9)=Q0(ii,1:2); %the orientation of the low b val eigVec1

            x0(3)=0.5;%  this is a case: just the center
            x0(4)=l1(ii);%D// lh_par
            x0(5)=mean([l2(ii) l3(ii)])/l1(ii);%Lh_par/Lh_per  --> Lh_per=x4*x5;
            x0(6)=2; %Dr_par of axon
            x0(7)=0.03; % this is a case from yaniv data
%              if check==1;
%        raw1=data(normal_b1,ii)/B0(ii); x0(3)=1;
%         Charmed_err(x0,raw1,Par)
%              end;
            [x1, resnorm] = lsqnonlin(@(x) Charmed_err(x,raw,Par),x0,lb,ub,options);

            x_1(:,ii)=x1;
            resnorm1(ii)=resnorm;
        end;

        res.x=x_1;
        res.resnorm=resnorm1;
        res.Par=Par;
        res.lb=lb;
        res.lb=ub;
        res.l0w_b_dif=A;
        res.brainMask=mask;
        res.brainInds=brainInds;
         save(outBaseName,'res')
      %  save(outBaseName,'x_1','resnorm1','A')
        fprintf('finish charmed fit and save it');

        return;

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

    case 'charmed_fitR'
        fprintf('Fitting %d Charmed model (SLOW!)...\n',nvox);

        %in order to get delta do that in a unix shell

        % dicomDumpHeader raw/dwi_g354_b10000_I0002.dcm
        % Stikov^Nikola  (M , 030Y): Date 20090519 Exam 1989, Series 5 , Image 2
        % 28 slices in series, Location: -6.411901474 (-8.4119  to 45.5881 )
        % (DWI CVs) b-value=10000, dwepi=354, diffusionGradDur=41124

        %[grads, maxB,delta,Delta] = dtiGradsBuildCharmed(maxG, Delta,delta,interact)
        %this is thecode we used to genarate the sequnse
        [grads, maxB,delta,Delta] = dtiGradsBuildCharmed([], [], delta,0);
        % maxG = 50.0; % mT/m * 1e-3/1e+6 = 1e-9
        % % For an optimized DTI sequence, we can get an effective G of srt(2)*maxG:
        % G = sqrt(2)*maxG;
        %bval = (2.*pi.*g).^2 * (G.*1e-9).^2 .* delta.^2 .* (Delta-delta./3)
        % bval = .8; % msec/micrometer^2 (8000 sec/mm^2 / 1000)
        % Delta = bval/((2*pi*g).^2 * (G*1e-9).^2 * delta.^2) + delta/3;
        %Delta=43.7;
        Del=(Delta-(delta/3));
        maxG = 50;
        G = grads*maxG;
        q =  G *1e-9 * g * delta;%  G*1e-9 -> T/m  g ->(cycles/millisecond)/T delta -> millisecond   q --> cycles
        % Check to make sure we computed q properly:
        %norm_q = sqrt(sum(q.^2,1))
        % b = (norm_q.*2*pi).^2*Del;
        % normB = sqrt(sum(b.^2,1))
        %

        Par(2)=4*pi^2;
        Par(3)=-Par(2)*Del; %need to check that this number is problematic!!!! what was it in the lowB Val?

        Par(1)=Delta;
        Par(5)=delta;
        Par(6)=Del;

        Dpar=1;
        Par(4)=Dpar;% micrometer^2/msec^2=1e-5 cm^2/sec^2  = micrometer^2/msec^2

        %     Par
        %     1- Delta
        %     2- 4Pi^2
        %     3- 4Pi^2*(Delta-(delta/3)
        %     4- 1=D_par
        %     5- delta
        %     6- (Delta-(delta/3));
        % q1=(bvecs*(G*1e-9)*g*delta*(2*pi))';

        %q1=(bvecs*(G*1e-9)*g*delta)';%*(2*pi))';
        %q units= T/um* cycles/millisecond)/T* ms*1/cycles =1/um
        %note  yaniv theta and phi are just the oposite of matlab so we take
        %them in reveres order. also matlab number running-  pi/2 <theta<pi/2
        %<pi<phi<phi in radians (for yanivs theta and phi- confusing!)
        %yaniv notition for Phi is the angel on the x-y serpase but he mesure
        %it from yaxes and matlab from x. so we must take it from pi/2
        %phi=pi/2-phi also for theta
        %from each phi number  and from theta .
        %q1=q/Delta;

        dwInds = bvals>0;
        B0 = mean(data(~dwInds,:),1);
        Dw = data(dwInds,:);

        %first use the compute one exponent dti with low bval:
        dt = dtiLoadDt6(fullfile(dt6_path,'dt6'));
        for dd=1:6
            dt6=squeeze(dt.dt6(:,:,:,dd));
            dtroi(dd,:)=dt6(ROICLmask);
        end;

        [eigVec,eigVal] = dtiEig(dtroi');
        [fa,md,rd] = dtiComputeFA(eigVal);
        l1 = eigVal(:,1);
        l2 = eigVal(:,2);
        l3 = eigVal(:,3);
        eig_L1(:,:)=eigVec(:,[1 2 3],1);
        %gof = zeros(1,nvox,'int16');

        lb=[.9  .1  .8  1e-6  1e-6]; %low bondery
        ub=[ 2   1   2  3     2]; %up bondery
        %  L1  L2/L1 L1ax a     b
        %Par(4)=Dpar*Del;

        options =  optimset('LevenbergMarquardt','on');

            % First compute the linear inversion matrix X^-1. This is
            % equivalent to the [-(1/b)] in the central DTI equation:
            % D = -(1/b) * log(S1/S0) (cf, Cercignani 2001; Basser 2002)
            % A = Xinv * log(data) : how we will represent it here
            %         Xinv = pinv(X);
            %         % To avoid log of zero warning:
            %         nz = data>0; logData = data; logData(nz) = log(data(nz));
            %         clear nz; % To avoid running out of memory.
            %         % Multiply Xinv * logData for each "page" or 2D matrix in a stack
            %         % (with numVoxels = number of "pages"). Each multiplication
            %         % involves Xinv (7 x nVols) * logData (nVols x 1)
            %         A = ndfun('mult', Xinv, logData); % Ainit = 7 x 1 x nVoxels
            % %         l=size(A);
            % %         A(l(1)+1,:,:)=1;
            %         normal_Tensor=A;
            %         [eigVec,eigVal] = dtiEig(squeeze(A(2:7,1,:))');
            %
            %         eig_L1(:,:)=eigVec(:,[1 2 3],1);
            %         l1 = eigVal(:,1);
            %         l2 = eigVal(:,2);
            %         l3 = eigVal(:,3);
            % Ainit = initial linear fit of seven
            % diffusion parameters for each voxel (6 directions + b0)
            %clear logData;
            %options = optimset('LargeScale','on','LevenbergMarquardt','on', 'Display', 'off', 'MaxIter', 50);
            %options = optimset('Display', 'off', 'MaxIter', 100);
            sigmaSq = sigma.^2;
            offset = 1e-6;

            %the matlab and charmed convention are opposing so we need to fix it.
            %first the use the opposit leter and then the mesure it in 90 degree
            %diference.
            [Q(:,1),Q(:,2) Q(:,3)] =cart2sph(q(1,:)',q(2,:)',q(3,:)');
            [Q0(:,1),Q0(:,2) Q0(:,3)] =cart2sph(eig_L1(:,1),eig_L1(:,2),eig_L1(:,3));%
            Q(:,1:2)=pi/2-Q(:,1:2);
            Q0(:,1:2)=pi/2-Q0(:,1:2);

            % we do this calculation ones as we assume that in the CC the angel is
            % fitted good.

            for j=1:length(Q0)

                Qper2(:,j)=Q(:,3).^2.*(1-(sin(Q(:,1)).*sin(Q0(j,1)).*cos(Q(:,2)-Q0(j,2))+cos(Q(:,1)).*cos(Q0(j,1))).^2);%Q+

                Qpar2(:,j)=Q(:,3).^2.*(sin(Q(:,1)).*sin(Q0(j,1)).*cos(Q(:,2)-Q0(j,2))+cos(Q(:,1)).*cos(Q0(j,1))).^2;
            end;
            %[Q1(:,2),Q1(:,1) Q1(:,3)] =cart2sph(q1(:,1),q1(:,2),q1(:,3));
            %     [THETA(:,1),PHI(:,1),R(:,1)] =cart2sph(q1(:,1),q1(:,2),q1(:,3));
            %
            %     Y_THETA(:,1)=pi/2-PHI(:,1);
            %     Y_PHI=(:,1)=pi/2-THETA(:,1);
            %
            % QRD2=Q(1,3)^2.*(1-(sin(Y_THETA(:,1)).*sin(x(1)).*cos(Y_PHI(:,2)-x(2))+cos(Y_THETA(:,1)).*cos(x(1))).^2);%Q+
            %
            % QLong2=Q(1,3)^2.*(sin(Y_THETA(:,1)).*sin(x(1)).*cos(Y_PHI(:,2)-x(2))+cos(Y_THETA(:,1)).*cos(x(1))).^2;%Q+
            %        Q(1,3)^2.*(sin(Q(:,1)).*sin(x(1)).*cos(Q(:,2)-x(2))+cos(Q(:,1)).*cos(x(1))).^2;%Q//

            % lb=ones(1,6)*offset ;
            % lb(3)=0.01;
            % lb(4)=0.01;
            % lb(5)=Par(4)*.1;
            % lb(6)=0.1;

            lb=[.9  .1  .8  1e-6  1e-6]; %low bondery
            ub=[ 2   1   2  3     2]; %up bondery
            %  L1  L2/L1 L1ax a     b
            %Par(4)=Dpar*Del;

            options =  optimset('LevenbergMarquardt','on');
            for(ii=1:nvox)
                raw=Dw(:,ii)/B0(ii);
                %[x0(2),x0(1)] =cart2sph(eig_L1(ii,1),eig_L1(ii,2),eig_L1(ii,3));%eigVecl1(vox(1),vox(2),vox(3),[1 2 3],1))%,....)
                %x0(1)=pi/2-x0(1);
                %x0(3)=pi/2-x0(3);
                x0(1)=l1(ii);%L// L_par
                x0(2)=mean([l2(ii) l3(ii)])/l1(ii);%L_par/L_per  --> Lper=x2*x1;
                x0(3)=2; %L_par of axon
                x0(4)=1; %alfa of gamma
                x0(5)=1; %beta of gamma
                %Par(1)=Fr(ii);

                if Wmcc(ii)>0;

                    %[x1, resnorm] = lsqnonlin(@(x) dtiRawCharmed_PD_BPD_Err(x,raw,Par,Qper2(dwInds,ii),Qpar2(dwInds,ii),Cmcc(ii),Wmcc(ii),Gmcc(ii),BPF(ii)),x0,lb,ub,options);
                    %[x1, resnorm] = lsqnonlin(@(x) dtiRawCharmed_PD_B_Err(x,raw,Par,Qper2(dwInds,ii),Qpar2(dwInds,ii),Cmcc(ii),BPF(ii)),x0,lb,ub,options);
                    [x1, resnorm] = lsqnonlin(@(x2) dtiRawCharmed_PD_Err(x2,raw,Dw(:,ii),B0(ii),Par,Qper2(dwInds,ii),Qpar2(dwInds,ii),Cmcc(ii),Wmcc(ii),Gmcc(ii)),x0,lb,ub,options);
                    Charmed_test_Education(x0,raw,Dw(:,ii),B0(ii),Par,Qper2(dwInds,ii),Qpar2(dwInds,ii),Cmcc(ii),Wmcc(ii),Gmcc(ii),Q(dwInds,:),Q0(ii,:),eig_L1(ii,:))
                    x_1(:,ii)=x1;
                    resnorm1(ii)=resnorm;
                end;
                %    [x(i,k,1:6),resnorm(i,k)] = lsqnonlin(@Charmedfun,x0,lb,ub,options) ;
            end;
            aMap=zeros(size(brainMask));
            bMap=zeros(size(brainMask));
            fitMap=zeros(size(brainMask));

            aMap(ROICLmask)=x_1(4,:);
            bMap(ROICLmask)=x_1(5,:);
            fitMap(ROICLmask)=resnorm1(:);

            save(outBaseName,'aMap','bMap','fitMap');
            % Par(1)=Fr(ii);
            return;
            % Par(2)=;4*pi^2;
            % Par(3)=-Par(2)*Del;
            %Par(4)=Dper;

            %x(1)=THETA
            % x(2)=PHI
            % x(3)=lpar
            % x(4)=lper
            % x(5)=Dpar
            % x(6)=R

        ;
            otherwise,
                error('unknown tensor fitting method "%s".',fitMethod);
        end

        if(bs.n>1)
            % the data use nvox*numVols doubles and each output uses nvox
            maxMemForStrides = bs.maxMem-(nvox*numVols+4*nvox)*8;
            % Xinv uses 7*numVols*bs.n doubles (8 bytes/dbl)
            % tmp uses 1*numVols*bs.n (=8*numVols*bs.n)
            % the resulting tensor fits use another 7*bs.n
            % We also allow for ~ 25% overhead per stride.
            stride = floor(maxMemForStrides./((8*numVols*bs.n+7*bs.n)*8*1.25));
            tic;
            fprintf('Running %d bootstrap samples on %d voxels (i.e. fitting %0.1f million tensors)- this may take a while!\n',bs.n,nvox,bs.n*nvox/1e6);
            % For the bootstrap, we set up a 3d Xinv matrix and then for each
            % voxel, solve for all the bootstrap tensor-fits at once using ndfun.
            if(isfield(bs,'permuteMatrix'))
                sampInds = zeros(numVols,bs.n);
                for(ii=1:numVols)
                    sampInds(ii,:) = bs.permuteMatrix{ii}(ceil(length(bs.permuteMatrix{ii}).*rand(1,bs.n)));
                end
            else
                % FIXME: allow for partial datasets, where numVols<nVolsPerRepeat*nRepeats.
                nRepeats = ceil(numVols/bs.nVolsPerRepeat);
                % The following will sample (with replacement) such that each repeat is
                % a coherent whole (ie. a full set of directions/bvals). You can check
                % this with mod(sampInds,nVolsPerRepeat).
                sampInds = repmat((1:bs.nVolsPerRepeat)',nRepeats,bs.n)+floor(nRepeats.*rand(bs.nVolsPerRepeat*nRepeats,bs.n))*bs.nVolsPerRepeat;
            end
            bsBvecs = bvecs(:,sampInds(:));
            bsBvals = bvals(sampInds(:));
            q = (bsBvecs.*sqrt(repmat(bsBvals./tau,3,1)))';
            sz = size(q);
            X = [ones(sz(1),1) -tau.*q(:,1).^2 -tau.*q(:,2).^2 -tau.*q(:,3).^2 -2*tau.*q(:,1).*q(:,2) -2*tau.*q(:,1).*q(:,3) -2*tau.*q(:,2).*q(:,3)];
            clear q bsBvals bsBvecs;
            X = permute(reshape(X',[7,numVols,bs.n]),[2 1 3]);
            %Xinv = ndfun('inv',X);
            fprintf('   Inverting %d X matrices...\n',bs.n);
            for(ii=1:size(X,3))
                Xinv(:,:,ii) = pinv(X(:,:,ii));
            end
            clear X;
            % To speed up the following loop, we do several (='stride') voxels per
            % iteration. If we had enough RAM, we could do them all in one go call
            % to ndfun to keep the voxel-loop in the c-code. But, the following
            % more pratical solution is a good compromise for a reasonably fast
            % solution (e.g., under 20 minutes for 500 bootstrap samples on ~200K
            % voxels.) Unfortunately, this optimization makes the following code
            % hard to read...
            Xinv = repmat(Xinv,[1,1,stride]);
            faStd = zeros(1,nvox);
            mdStd = zeros(1,nvox);
            pddDisp = zeros(1,nvox);
            fprintf('   Multiplying log(DWI) by Xinv (%d voxel stride)...\n',stride);
            for(ii=1:stride:nvox)
                sInd = ii; eInd = min(ii+stride-1,nvox);
                n = eInd-sInd+1;
                tmp = reshape(data(sampInds(:),1,sInd:eInd),[numVols,1,bs.n*n]);
                t = ndfun('mult', Xinv(:,:,1:size(tmp,2)), tmp);
                t = squeeze(t(2:7,:,:));
                t = reshape(t,6,bs.n,n);
                % We now have bs.n tensors- use them to estimate tensor fit
                % variances like FA, MD and PDD variance.
                %
                % dtiEig likes the 6 tensor vals to be in the 2nd (or 4th) dim.
                % we'll permute so that the n voxels for this stride are in the
                % first dim and the bs.n are in the 3rd dim.
                %
                % *** WORK HERE- insert better PDD dispersion fitting here (e.g.
                % Bingham distribution)
                [vec,val] = dtiEig(permute(t,[3 1 2]));
                % Convert vec to an Mx3xN array of PDD vectors where N is the
                % bootstrap dim- the dim along which we'll collapse to compute
                % dispersions.
                vec = permute(vec(:,:,1,:),[1 2 4 3]);
                %badVals = any(val<0,2);
                [fa,md] = dtiComputeFA(val);
                faStd(sInd:eInd) = std(fa,0,2);
                %         keyboard;
                %         tic
                %         for(jj=1:size(fa,1))
                %             [pdfY,pdfX] = ksdensity(fa(jj,:),'function','cdf');
                %             xi = interp1(pdfY+rand(size(pdfY))*1e-9,pdfX,[0.025 0.975]);
                %             faLb(sInd+jj-1) = xi(1); faUb(sInd+jj-1) = xi(2);
                %         end
                %         toc
                mdStd(sInd:eInd) = std(md,0,2);
                clear t val;
                [junk,pddDisp(sInd:eInd)] = dtiDirMean(vec);
            end
            clear Xinv tmp vec;
            % Convert dispersion to angle in degrees
            % We get a few values just below zero in edge voxels, so we clip.
            pddDisp(pddDisp<0) = 0;
            pddDisp = asin(sqrt(pddDisp))./pi.*180;
            fprintf('   Elapsed time: %f minutes.\n',toc/60);
            tmp = zeros(size(brainMask));
            tmp(brainInds) = faStd; faStd = tmp;
            tmp(brainInds) = mdStd; mdStd = tmp;
            tmp(brainInds) = pddDisp; pddDisp = tmp;
        else
            faStd = [];
            mdStd = [];
            pddDisp = [];
        end

        tn = tempname;
        disp(['saving temp data to ' tn '...']);
        save(tn);

        % Computing the mean b0 (in log space) is exactly equivalent to pulling the
        % b0 from the model fit (A(1,:,:)). We use the log-mean b0 that we already
        % computed since it saves us a little computation (one less exp) and it has
        % no background voxels masked away.
        %b0 = zeros(size(brainMask));
        %b0(brainInds) = exp(squeeze(A(1,:,:)));
        tmp = zeros(size(brainMask));
        dt6 = zeros([size(brainMask),6]);
        for(ii=1:6)
            tmp(brainInds) = squeeze(A(ii+1,:,:));
            dt6(:,:,:,ii) = tmp;
        end

        %% Check ADC units and adjust if necessary
        %
        if(isempty(adcUnits))
            % Always convert to our standard ADC units (micrometer^2/msec)
            [curUnitStr,scale,adcUnits] = dtiGuessDiffusivityUnits(dt6);
            if(scale~=1)
                fprintf('Converting %s to %s with scale = %f.\n',curUnitStr,adcUnits,scale);
                dt6 = dt6.*scale;
            end
        end

        %% Compute a rough white-matter mask
        %
        [fa,md] = dtiComputeFA(dt6);
        wmMask = brainMask & fa>.15 & (md<1.1 | fa>0.4);
        wmMask = dtiCleanImageMask(wmMask,0,0);

        %% Save all results
        %
        if(~exist(outBaseName,'dir'))
            mkdir(outBaseName);
        end
        if(~exist(binDirName,'dir'))
            mkdir(binDirName);
        end
        params.nBootSamps = bs.n;
        params.buildDate = datestr(now,'yyyy-mm-dd HH:MM');
        l = license('inuse');
        params.buildId = sprintf('%s on Matlab R%s (%s)',l(1).user,version('-release'),computer);
        params.rawDataDir = dataDir;
        % We assume that the raw data file is a directory inside the 'subject'
        % directory.
        params.subDir = fileparts(dataDir);

        % We want all the important file names to be relative so that they are
        % platform-independent. The only platform-dependent path should be
        % 'homeDir'. As long as the dt6 project file stays in the same dir as the
        % bin dir, we shouldn't need 'homeDir' to find everything.
        %
        %  TONY: Every filename now is relative to the directory above the
        %  directory containing this dt6.  Thus we have no system dependent
        %  information.  Just make sure the relative paths stay the same.
        %[files.homeDir,files.binDir] = fileparts(binDirName);
        [fullParentDir, binDir] = fileparts(binDirName);
        [ppBinDir, pBinDir] = fileparts(fullParentDir);
        pBinDir = fullfile(pBinDir,binDir);
        files.b0 = fullfile(pBinDir,'b0.nii.gz');
        files.brainMask = fullfile(pBinDir,'brainMask.nii.gz');
        files.wmMask = fullfile(pBinDir,'wmMask.nii.gz');
        files.tensors = fullfile(pBinDir,'tensors.nii.gz');
        % description can have up to 80 chars
        desc = [params.buildDate ' ' params.buildId];
        if(length(desc)>80), disp('NOTE: description field clipped to 80 chars.'); end
        dtiWriteNiftiWrapper(int16(round(b0)), xformToAcPc, fullfile(ppBinDir,files.b0), 1, desc, 'b0');
        dtiWriteNiftiWrapper(uint8(brainMask), xformToAcPc, fullfile(ppBinDir,files.brainMask), 1, desc, 'brainMask');
        dtiWriteNiftiWrapper(uint8(wmMask), xformToAcPc, fullfile(ppBinDir,files.wmMask), 1, desc, 'whiteMatterMask');
        if(~isempty(gof))
            tmp = zeros(size(brainMask),'int16');
            tmp(brainInds) = gof;
            gof = tmp;
            tmpVol = zeros([size(brainMask),numVols],'uint8');
            tmp = zeros(size(brainMask),'uint8');
            for(ii=1:numVols)
                tmp(brainInds) = outliers(ii,:);
                tmpVol(:,:,:,ii) = tmp;
            end
            outliers = tmpVol;
            files.gof = fullfile(pBinDir,'gof.nii.gz');
            files.outliers = fullfile(pBinDir,'outliers.nii.gz');
            dtiWriteNiftiWrapper(gof, xformToAcPc, fullfile(ppBinDir,files.gof), 1, desc, 'GOF');
            dtiWriteNiftiWrapper(outliers, xformToAcPc, fullfile(ppBinDir,files.outliers), 1, desc, 'outlier mask');
            %Create summary image of outliers.nii.gz that can be viewed as an image
            % when loaded into DTIfiberUI.
            %outlierImage=niftiRead(fullfile(pBinDir,files.outliers));
            outlierImage=niftiRead(files.outliers);
            outlierImage.data=sum(outlierImage.data,4);
            outlierImage.fname = fullfile(ppBinDir,pBinDir,'outlier_sum_image.nii.gz');
            writeFileNifti(outlierImage);
        end
        % NIFTI convention is for the 6 unique tensor elements stored in the 5th
        % dim in lower-triangular, row-order (Dxx Dxy Dyy Dxz Dyz Dzz). NIFTI
        % reserves the 4th dim for time, so in the case of a time-invatiant tensor,
        % we just leave a singleton 4th dim. Our own internal convention is
        % [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz], so we use the code below to convert to
        % the NIFTI order and dt6=squeeze(ni.data(:,:,:,1,[1 3 6 2 4 5])); to get
        % back to our convention. FOr reference- the 3x3 tensor matrix is:
        %    Dxx Dxy Dxz
        %    Dxy Dyy Dyz
        %    Dxz Dyz Dzz
        dt6 = dt6(:,:,:,[1 4 2 5 6 3]);
        sz = size(dt6);
        dt6 = reshape(dt6,[sz(1:3),1,sz(4)]);
        dtiWriteNiftiWrapper(dt6, xformToAcPc, fullfile(ppBinDir,files.tensors), 1, desc, ['DTI ' adcUnits]);
        if(bs.n>1)
            files.faStd = fullfile(pBinDir,'faStd.nii.gz');
            files.mdStd = fullfile(pBinDir,'mdStd.nii.gz');
            files.pddDisp = fullfile(pBinDir,'pddDispersion.nii.gz');
            dtiWriteNiftiWrapper(single(faStd), xformToAcPc, fullfile(ppBinDir,files.faStd), 1, desc, 'FA stdev');
            dtiWriteNiftiWrapper(single(mdStd), xformToAcPc, fullfile(ppBinDir,files.mdStd), 1, desc, 'MD stdev');
            dtiWriteNiftiWrapper(pddDisp, xformToAcPc, fullfile(ppBinDir,files.pddDisp), 1, desc, 'PDD disp (deg)');
        end
        save(dt6FileName,'adcUnits','params','files');
        disp('Finished writing dt6 file--line634 of dtiRawFitTensor');
        if(nargout<1), clear dt6; end

        return;

        % % To run this on a bunch of subjects
        % bd = '/biac3/wandell4/data/reading_longitude/dti_y1234';
        % rf = 'dti_g13_b800_aligned.';
        % of = 'dti06rt';
        % d = dir(fullfile(bd,'*0*'));
        % for(ii=1:numel(d))
        %     sd = fullfile(bd,d(ii).name);
        %     dwRaw = fullfile(sd,'raw',[rf 'nii.gz']);
        %     dwBvc = fullfile(sd,'raw',[rf 'bvecs']);
        %     dwBvl = fullfile(sd,'raw',[rf 'bvals']);
        %     out = fullfile(sd,of);
        %     if(exist(sd,'dir') && exist(dwRaw,'file') && exist(dwBvc,'file') && exist(dwBvl,'file') && ~exist(out,'dir'))
        %         fprintf('Processing %s (%d of %d)...\n',sd,ii,numel(d));
        %         dtiRawFitTensor(dwRaw, dwBvc, dwBvl, out, [], 'rt');
        %     end
        % end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function s = SetDefaultParams(s);
% s = SetDefaultParams(s);
% Sets default parameters
% s: user-specified parameters that are used instead of defaults

if (isfield(s, 'gamma') == 0),
    s.gamma = 1;
end

if (isfield(s, 'beta') == 0),
    s.beta = 1;
end

if (isfield(s, 'const_factor') == 0),
    s.const_factor = 40;
end

if (isfield(s, 'type4_rank') == 0),
    s.type4_rank = 5;
end

if (isfield(s, 'thresh') == 0),
   s.thresh = 10e-3;
end

if (isfield(s, 'k') == 0),
    s.k = 4;
end

if (isfield(s, 'max_iters') == 0),
    s.max_iters = 100000;
end


text:
function varargout = subdir(varargin)
%SUBDIR Performs a recursive file search
%
% subdir
% subdir(name)
% files = subdir(...)
%
% This function performs a recursive file search.  The input and output
% format is identical to the dir function.
%
% Input variables:
%
%   name:   pathname or filename for search, can be absolute or relative
%           and wildcards (*) are allowed.  If ommitted, the files in the
%           current working directory and its child folders are returned
%
% Output variables:
%
%   files:  m x 1 structure with the following fields:
%           name:   full filename
%           date:   modification date timestamp
%           bytes:  number of bytes allocated to the file
%           isdir:  1 if name is a directory; 0 if no
%
% Example:
%
%   >> a = subdir(fullfile(matlabroot, 'toolbox', 'matlab', '*.mat'))
%
%   a =
%
%   67x1 struct array with fields:
%       name
%       date
%       bytes
%       isdir
%
%   >> a(2)
%
%   ans =
%
%        name: '/Applications/MATLAB73/toolbox/matlab/audiovideo/chirp.mat'
%        date: '14-Mar-2004 07:31:48'
%       bytes: 25276
%       isdir: 0
%
% See also:
%
%   dir

% Copyright 2006 Kelly Kearney

%---------------------------
% Get folder and filter
%---------------------------

narginchk(0,1);
nargoutchk(0,1);

if nargin == 0
    folder = pwd;
    filter = '*';
else
    [folder, name, ext] = fileparts(varargin{1});
    if isempty(folder)
        folder = pwd;
    end
    if isempty(ext)
        if isdir(fullfile(folder, name))
            folder = fullfile(folder, name);
            filter = '*';
        else
            filter = [name ext];
        end
    else
        filter = [name ext];
    end
    if ~isdir(folder)
        error('Folder (%s) not found', folder);
    end
end

%---------------------------
% Search all folders
%---------------------------

pathstr = genpath_local(folder);
pathfolders = regexp(pathstr, pathsep, 'split');  % Same as strsplit without the error checking
pathfolders = pathfolders(~cellfun('isempty', pathfolders));  % Remove any empty cells

Files = [];
pathandfilt = fullfile(pathfolders, filter);
for ifolder = 1:length(pathandfilt)
    NewFiles = dir(pathandfilt{ifolder});
    if ~isempty(NewFiles)
        fullnames = cellfun(@(a) fullfile(pathfolders{ifolder}, a), {NewFiles.name}, 'UniformOutput', false);
        [NewFiles.name] = deal(fullnames{:});
        Files = [Files; NewFiles];
    end
end

%---------------------------
% Prune . and ..
%---------------------------

if ~isempty(Files)
    [~, ~, tail] = cellfun(@fileparts, {Files(:).name}, 'UniformOutput', false);
    dottest = cellfun(@(x) isempty(regexp(x, '\.+(\w+$)', 'once')), tail);
    Files(dottest & [Files(:).isdir]) = [];
end

%---------------------------
% Output
%---------------------------

if nargout == 0
    if ~isempty(Files)
        fprintf('\n');
        fprintf('%s\n', Files.name);
        fprintf('\n');
    end
elseif nargout == 1
    varargout{1} = Files;
end

function [p] = genpath_local(d)
% Modified genpath that doesn't ignore:
%     - Folders named 'private'
%     - MATLAB class folders (folder name starts with '@')
%     - MATLAB package folders (folder name starts with '+')

files = dir(d);
if isempty(files)
  return
end
p = '';  % Initialize output

% Add d to the path even if it is empty.
p = [p d pathsep];

% Set logical vector for subdirectory entries in d
isdir = logical(cat(1,files.isdir));
dirs = files(isdir);  % Select only directory entries from the current listing

for i=1:length(dirs)
   dirname = dirs(i).name;
   if    ~strcmp( dirname,'.') && ~strcmp( dirname,'..')
       p = [p genpath(fullfile(d,dirname))];  % Recursive calling of this function.
   end
end

repo: CaImAn-MATLAB