text:
function [M_final,shifts,template,options,col_shift] = normcorre(Y,options,template)

% online motion correction through DFT subpixel registration
% Based on the dftregistration.m function from Manuel Guizar and Jim Fienup

% INPUTS
% Y:                Input data, can be already loaded in memory as a 3D
%                   tensor, a memory mapped file, or a pointer to a tiff stack
% options:          options structure for motion correction (optional, rigid registration is performed if not provided)
% template:         provide template (optional)

% OUTPUTS
% M_final:          motion corrected data
% shifts:           originally calculated shifts
% template:         calculated template
% options:          options structure (if modified)
% col_shift:        relative shift due to bi-directional scanning

%% first determine filetype

nd = 2 + (options.d3 > 1); %max(length(sizY)-1,2);                    % determine whether imaging is 2d or 3d

if isa(Y,'char')
    [~,~,ext] = fileparts(Y);
    ext = ext(2:end);
    if strcmpi(ext,'tif') || strcmpi(ext,'tiff')
        tiffInfo = imfinfo(Y);
        filetype = 'tif';
        T = length(tiffInfo);
        if nd == 3
            sizY = [tiffInfo(1).Height,tiffInfo(1).Width,T,1];
        else
            sizY = [tiffInfo(1).Height,tiffInfo(1).Width,T];
        end
    elseif strcmpi(ext,'mat')
        filetype = 'mem';
        Y = matfile(Y,'Writable',true);
        details = whos(Y);
        var_sizes = [details.bytes];
        [~,var_ind] = max(var_sizes);
        var_name = details(var_ind).name;
        sizY = size(Y,var_name);
    elseif strcmpi(ext,'hdf5') || strcmpi(ext,'h5')
        filetype = 'hdf5';
        fileinfo = hdf5info(Y);
        sizY = fileinfo.GroupHierarchy.Datasets.Dims;
    end
elseif isobject(Y)
    filetype = 'mem';
    var_name = 'Y';
    sizY = size(Y,var_name);
else % array loaded in memory
    filetype = 'mat';
    Y = single(Y);
    sizY = size(Y);
end

if length(sizY) == nd
    T = 1;
else
    T = sizY(nd+1);
end
sizY = sizY(1:nd);

%% set default parameters if not present

if ~exist('options','var') || isempty(options)
    options = NoRMCorreSetParms('d1',sizY(1),'d2',sizY(2));
    if nd > 2; options.d3 = sizY(3); end
end

options.bin_width = min(options.bin_width,T+1);
options.mem_batch_size = min(options.mem_batch_size,T);
options.buffer_width = min(options.buffer_width,ceil(T/options.bin_width));
options.init_batch = min(options.init_batch,T);

memmap = options.memmap;
grid_size = options.grid_size;
mot_uf = options.mot_uf;
min_patch_size = options.min_patch_size;
overlap_pre = options.overlap_pre;
overlap_post = options.overlap_post;
upd_template = options.upd_template;
bin_width = options.bin_width;
buffer_width = options.buffer_width;
max_dev = options.max_dev;
init_batch = options.init_batch;
us_fac = options.us_fac;
method = options.method;
plot_flag = options.plot_flag*(nd==2);
filename = options.mem_filename;
use_parallel = options.use_parallel;
make_avi = options.make_avi;
name = options.name;
fr = options.fr;
iter = options.iter;
add_value = options.add_value;
max_shift = options.max_shift;
print_msg = options.print_msg;
if strcmpi(options.boundary,'nan')
    fill_value = NaN;
else
    fill_value = add_value;
end

%% first check for offset due to bi-directional scanning

if options.correct_bidir && isempty(options.col_shift)
    col_shift = correct_bidirectional_offset(Y,options.nFrames,options.bidir_us);
elseif ~isempty(options.col_shift)
    col_shift = options.col_shift;
else
    col_shift = 0;
end
options.col_shift = col_shift;
if col_shift
    if print_msg; fprintf('Offset %1.1d pixels due to bidirectional scanning detected. \n',col_shift); end
    if strcmpi(options.shifts_method,'fft')
        options.shifts_method = 'cubic';
        if print_msg; fprintf('Cubic shifts will be applied. \n'); end
    end
end
%% read initial batch and compute template
init_batch = min(T,init_batch);
interval = ceil(T/2-init_batch/2+1):floor(T/2+init_batch/2);
if exist('template','var')
    init_batch = min(init_batch,1);
end
switch filetype
    case 'tif'
        Y_temp = read_file(Y,interval(1),init_batch,[],tiffInfo);
    case 'hdf5'
        Y_temp = read_file(Y,interval(1),init_batch);
    case 'mem'
        if nd == 2; Y_temp = Y.(var_name)(:,:,interval); elseif nd == 3; Y_temp = Y.(var_name)(:,:,:,interval); end
    case 'mat'
        if nd == 2; Y_temp = Y(:,:,interval); elseif nd == 3; Y_temp = Y(:,:,:,interval); end
end

data_type = class(Y_temp);
Y_temp = single(Y_temp);

if nargin < 3 || isempty(template)
    if print_msg; fprintf('Registering the first %i frames just to obtain a good template....',init_batch); end
    template_in = median(Y_temp,nd+1)+add_value;
    fftTemp = fftn(template_in);
    for t = 1:size(Y_temp,nd+1)
        if nd == 2
            [~,Greg] = dftregistration_min_max(fftTemp,fftn(Y_temp(:,:,t)),us_fac,-max_shift,max_shift,options.phase_flag);
        end
        if nd == 3
            [~,Greg] = dftregistration_min_max_3d(fftTemp,fftn(Y_temp(:,:,:,t)),us_fac,-max_shift,max_shift,options.phase_flag);
        end
        M_temp = real(ifftn(Greg));
        template_in = template_in*(t-1)/t + M_temp/t;
    end
    template_in = template_in + add_value;
    if print_msg; fprintf('..done. \n'); end
else
    template_in = single(template + add_value);
end

[d1,d2,d3,~] = size(Y_temp);
if nd == 2; d3 = 1; end
%% setup grids for patches

[xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,xx_us,xx_uf,yy_us,yy_uf,zz_us,zz_uf] = construct_grid(grid_size,mot_uf,d1,d2,d3,min_patch_size);
shifts = struct('shifts',cell(T,1),'shifts_up',cell(T,1),'diff',cell(T,1));
temp_cell = mat2cell_ov(template_in,xx_us,xx_uf,yy_us,yy_uf,zz_us,zz_uf,overlap_post,sizY);

%% precompute some quantities that are used repetitively for template matching and applying shifts
Nr = cell(size(temp_cell));
Nc = cell(size(temp_cell));
Np = cell(size(temp_cell));
Bs = cell(size(temp_cell));
for i = 1:length(xx_us)
    for j = 1:length(yy_us)
        for k = 1:length(zz_us)
            [nr,nc,np] = size(temp_cell{i,j,k});
            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);
            [Nc{i,j,k},Nr{i,j,k},Np{i,j,k}] = meshgrid(nc,nr,np);
            extended_grid = [max(xx_us(i)-overlap_post(1),1),min(xx_uf(i)+overlap_post(1),d1),max(yy_us(j)-overlap_post(2),1),min(yy_uf(j)+overlap_post(2),d2),max(zz_us(k)-overlap_post(3),1),min(zz_uf(k)+overlap_post(3),d3)];
            Bs{i,j,k} = permute(construct_weights([xx_us(i),xx_uf(i),yy_us(j),yy_uf(j),zz_us(k),zz_uf(k)],extended_grid),[2,1,3]);
        end
    end
end
if nd == 2; Np = cellfun(@(x) 0,Nr,'un',0); end

%%
%maxNumCompThreads(2);
template = mat2cell_ov(template_in,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY);
temp_mat = template_in;
fftTemp = cellfun(@fftn,template,'un',0);
fftTempMat = fftn(temp_mat);
if nd == 2; buffer = mat2cell_ov(zeros(d1,d2,bin_width,'single'),xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end
if nd == 3; buffer = mat2cell_ov(zeros(d1,d2,d3,bin_width,'single'),xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end

if ~strcmpi(options.output_type,'mat')
    %options.mem_batch_size = min(round(options.mem_batch_size/bin_width)*bin_width,T);
    if nd == 2; mem_buffer = zeros(d1,d2,options.mem_batch_size,'single'); end
    if nd == 3; mem_buffer = zeros(d1,d2,d3,options.mem_batch_size,'single'); end
end

switch lower(options.output_type)
    case 'mat'
        M_final = zeros([sizY,T],data_type);
    case 'memmap'
        M_final = matfile(filename,'Writable',true);
        if nd == 2; M_final.Y(d1,d2,T) = zeros(1,data_type); end
        if nd == 3; M_final.Y(d1,d2,d3,T) = zeros(1,data_type); end
        M_final.Yr(d1*d2*d3,T) = zeros(1,data_type);
    case {'hdf5','h5'}
        if exist(options.h5_filename,'file')
            [pathstr,fname,ext] = fileparts(options.h5_filename);
            new_filename = fullfile(pathstr,[fname,'_',datestr(now,30),ext]);
            warning_msg = ['File ',options.h5_filename,'already exists. Saving motion corrected file as',new_filename];
            warning('%s',warning_msg);
            options.h5_filename = new_filename;
        end
        M_final = options.h5_filename;
        if nd == 2
            h5create(options.h5_filename,['/',options.h5_groupname],[d1,d2,Inf],'Chunksize',[d1,d2,options.mem_batch_size],'Datatype',data_type);
        elseif nd == 3
            h5create(options.h5_filename,['/',options.h5_groupname],[d1,d2,d3,Inf],'Chunksize',[d1,d2,d3,options.mem_batch_size],'Datatype',data_type);
        end
    case {'tif','tiff'}
        M_final = options.tiff_filename;
        opts_tiff.append = true;
        opts_tiff.big = true;
        if nd == 3
            error('Saving volumetric tiff stacks is currently not supported. Use a different filetype');
        end
    otherwise
        error('This filetype is currently not supported')
end

%%
if plot_flag
    if make_avi
        vidObj = VideoWriter(name);
        set(vidObj,'FrameRate',fr);
        open(vidObj);
    end
    if strcmpi(filetype,'mat')
        nnY = quantile(Y(:),0.005);
        mmY = quantile(Y(:),0.995);
    else
        nnY = quantile(Y_temp(:),0.005);
        mmY = quantile(Y_temp(:),0.995);
    end
    fig = figure;
        screensize = get(0,'Screensize' );
        fac = min(min((screensize(3:4)-200)./[d2,d1]),10);
        set(gcf, 'PaperUnits', 'points', 'Units', 'points');
        set(gcf, 'Position', round([100 100 fac*d2 fac*d1]));
end
cnt_buf = 0;
if print_msg; fprintf('Template initialization complete.  Now registering all the frames with new template. \n'); end
%%
prevstr = [];
for it = 1:iter
    if it < iter; plot_flag = 0; else plot_flag = options.plot_flag; end
    for t = 1:T
        switch filetype
            case 'tif'
                Yt = single(imread(Y,'Index',t,'Info',tiffInfo));
            case 'hdf5'
                Yt = single(h5read(Y,'/mov',[ones(1,nd),t],[sizY(1:nd),1]));
            case 'mem'
                if nd == 2; Yt = single(Y.(var_name)(:,:,t)); end
                if nd == 3; Yt = single(Y.(var_name)(:,:,:,t)); end
            case 'mat'
                if nd == 2; Yt = single(Y(:,:,t)); end
                if nd == 3; Yt = single(Y(:,:,:,t)); end
        end
        minY = min(Yt(:));
        maxY = max(Yt(:));
        Yt = Yt + add_value;
        ind = rem(t,bin_width) + bin_width*(rem(t,bin_width)==0);
        Yc = mat2cell_ov(Yt,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY);
        fftY = cellfun(@fftn, Yc, 'un',0);

        M_fin = cell(length(xx_us),length(yy_us),length(zz_us)); %zeros(size(Y_temp));
        shifts_temp = zeros(length(xx_s),length(yy_s),length(zz_s),nd);
        diff_temp = zeros(length(xx_s),length(yy_s),length(zz_s));
        if numel(M_fin) > 1
            if nd == 2; out_rig = dftregistration_min_max(fftTempMat,fftn(Yt),us_fac,-max_shift,max_shift,options.phase_flag); lb = out_rig(3:4); ub = out_rig(3:4); end
            if nd == 3; out_rig = dftregistration_min_max_3d(fftTempMat,fftn(Yt),1,-max_shift,max_shift,options.phase_flag); lb = out_rig(3:5); ub = out_rig(3:5); end
        else
            lb = -max_shift(1,nd);
            ub = max_shift(1,nd);
            max_dev = 0*max_dev;
        end
        if ~use_parallel
            for i = 1:length(xx_s)
                for j = 1:length(yy_s)
                    for k = 1:length(zz_s)
                        if nd == 2
                            %[output,Greg] = dftregistration_max(fftTemp{i,j,k},fftY{i,j,k},us_fac,max_shift);
                            [output,Greg] = dftregistration_min_max(fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev(1:2),ub+max_dev(1:2),options.phase_flag);
                        elseif nd == 3
                            %[output,Greg] = dftregistration_max_3d(fftTemp{i,j,k},fftY{i,j,k},us_fac,max_shift);
                            [output,Greg] = dftregistration_min_max_3d(fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev,ub+max_dev,options.phase_flag);
                            shifts_temp(i,j,k,3) = output(5);
                        end
                        M_temp = real(ifftn(Greg));
                        M_temp = remove_boundaries(M_temp,output(3:end),'copy',template{i,j,k});
                        if nd == 2; buffer{i,j,k}(:,:,ind) = M_temp;  end
                        if nd == 3; buffer{i,j,k}(:,:,:,ind) = M_temp;  end
                        shifts_temp(i,j,k,1) = output(3);
                        shifts_temp(i,j,k,2) = output(4);
                        diff_temp(i,j,k) = output(2);
                        if all([length(xx_s),length(yy_s),length(zz_s)] == 1)
                            M_fin{i,j,k} = remove_boundaries(M_temp,output(3:end),options.boundary,template{i,j,k},add_value);
                        end
                    end
                end
            end
        else
            Mt2 = cell(length(xx_s)*length(yy_s)*length(zz_s),1);
            shifts_cell = cell(length(xx_s)*length(yy_s)*length(zz_s),1);
            diff_cell = cell(length(xx_s)*length(yy_s)*length(zz_s),1);
            for ii = length(xx_s)*length(yy_s)*length(zz_s):-1:1
                [i,j,k] = ind2sub([length(xx_s),length(yy_s),length(zz_s)],ii);
                if nd == 2; future_results(ii) = parfeval(@dftregistration_min_max,2,fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev(1:2),ub+max_dev(1:2),options.phase_flag); end
                if nd == 3; future_results(ii) = parfeval(@dftregistration_min_max_3d,2,fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev,ub+max_dev,options.phase_flag); end
            end
            for i = 1:length(xx_s)*length(yy_s)*length(zz_s)
                [ii,output,Greg] = fetchNext(future_results);
                M_temp = real(ifftn(Greg));
                Mt2{ii} = M_temp;
                shifts_cell{ii} = output(3:end);
                diff_cell{ii} = output(2);
            end
%             parfor ii = 1:length(xx_s)*length(yy_s)*length(zz_s)
%                 [i,j,k] = ind2sub([length(xx_s),length(yy_s),length(zz_s)],ii);
%                 %if nd == 2; [output,Greg] = dftregistration_max(fftTemp{i,j,k},fftY{i,j,k},us_fac,max_shift); end
%                 %if nd == 3; [output,Greg] = dftregistration_max_3d(fftTemp{i,j,k},fftY{i,j,k},us_fac,max_shift); end
%                 if nd == 2; [output,Greg] = dftregistration_min_max(fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev(1:2),ub+max_dev(1:2),options.phase_flag); end
%                 if nd == 3; [output,Greg] = dftregistration_min_max_3d(fftTemp{i,j,k},fftY{i,j,k},us_fac,lb-max_dev,ub+max_dev,options.phase_flag); end
%                 M_temp = real(ifftn(Greg));
%                 Mt2{ii} = M_temp;
%                 shifts_cell{ii} = output(3:end);
%                 diff_cell{ii} = output(2);
%             end
            for ii = 1:length(xx_s)*length(yy_s)*length(zz_s)
                 [i,j,k] = ind2sub([length(xx_s),length(yy_s),length(zz_s)],ii);
                 if nd == 2; buffer{i,j,k}(:,:,ind) = Mt2{ii}; end
                 if nd == 3; buffer{i,j,k}(:,:,:,ind) = Mt2{ii}; end
                 if mot_uf == 1
                     M_fin{i,j,k} = Mt2{ii};
                 end
                 shifts_temp(i,j,k,:) = shifts_cell{ii};
                 diff_temp(i,j,k) = diff_cell{ii};
            end
        end
        shifts(t).shifts = shifts_temp;
        shifts(t).diff = diff_temp;

        switch lower(options.shifts_method)
            case 'fft'
                if any([length(xx_s),length(yy_s),length(zz_s)] > 1)
                    if ~isfield(options,'shifts_method'); options.shifts_method = 'FFT'; end
                    if mot_uf(3) > 1
                        do = [length(xx_us),length(yy_us),length(zz_us)]./[length(xx_s),length(yy_s),length(zz_s)];
                        ds = [length(xx_s),length(yy_s),length(zz_s)];
                        dim = [length(xx_us),length(yy_us),length(zz_us)];
                        [Xq,Yq,Zq] = meshgrid(linspace((1+1/do(2))/2,ds(2)+(1-1/do(2))/2,dim(2)),linspace((1+1/do(1))/2,ds(1)+(1-1/do(1))/2,dim(1)),linspace((1+1/do(3))/2,ds(3)+(1-1/do(3))/2,dim(3)));
                        %tform = affine3d(diag([mot_uf(:);1]));
                        %tform = affine3d(diag([mot_uf([2,1,3])';1]));
                        %diff_up = imwarp(diff_temp,tform,'OutputView',imref3d([length(xx_uf),length(yy_uf),length(zz_uf)]));
                        %diff_up = imwarp(diff_temp,tform,'OutputView',imref3d([length(xx_uf),length(yy_uf),length(zz_uf)]),'SmoothEdges',true);
                        diff_up = interp3(diff_temp,Xq,Yq,Zq,'makima');
                        shifts_up = zeros([size(diff_up),3]);
                        %for dm = 1:3; shifts_up(:,:,:,dm) = imwarp(shifts_temp(:,:,:,dm),tform,'OutputView',imref3d([length(xx_uf),length(yy_uf),length(zz_uf)]),'SmoothEdges',true); end
                        for dm = 1:3; shifts_up(:,:,:,dm) = interp3(shifts_temp(:,:,:,dm),Xq,Yq,Zq,'makima'); end
                    else
                        shifts_up = imresize(shifts_temp,[length(xx_uf),length(yy_uf)]);
                        diff_up = imresize(diff_temp,[length(xx_uf),length(yy_uf)]);
                    end

                    shifts(t).shifts_up = shifts_up;
                    shifts(t).diff = diff_up;
                    for i = 1:length(xx_uf)
                        for j = 1:length(yy_uf)
                            for k = 1:length(zz_uf)
                                extended_grid = [max(xx_us(i)-overlap_post(1),1),min(xx_uf(i)+overlap_post(1),d1),max(yy_us(j)-overlap_post(2),1),min(yy_uf(j)+overlap_post(2),d2),max(zz_us(k)-overlap_post(3),1),min(zz_uf(k)+overlap_post(3),d3)];
                                I_temp = Yt(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6));
                                M_fin{i,j,k} = shift_reconstruct(I_temp,shifts_up(i,j,k,:),diff_up(i,j,k),us_fac,Nr{i,j,k},Nc{i,j,k},Np{i,j,k},options.boundary,add_value);
                                %M_fin{i,j,k} = shift_reconstruct2(I_temp,shifts_up(i,j,k,:),'bilinear',diff_up(i,j,k),us_fac,Nr{i,j,k},Nc{i,j,k},Np{i,j,k},options.boundary,add_value);
                            end
                        end
                    end
                else
                    shifts_up = shifts_temp;
                    shifts(t).shifts_up = shifts(t).shifts;
                end
                gx = max(abs(reshape(diff(shifts_up,[],1),[],1)));
                gy = max(abs(reshape(diff(shifts_up,[],2),[],1)));
                gz = max(abs(reshape(diff(shifts_up,[],3),[],1)));
                flag_interp = max([gx;gy;gz;0])<0.5;      % detect possible smearing

                if flag_interp
                    Mf = cell2mat_ov_sum(M_fin,xx_us,xx_uf,yy_us,yy_uf,zz_us,zz_uf,overlap_post,sizY,Bs) - add_value;
                else
                    Mf = cell2mat_ov(M_fin,xx_us,xx_uf,yy_us,yy_uf,zz_us,zz_uf,overlap_post,sizY) - add_value;
                end

            otherwise
                shifts(t).shifts_up = shifts(t).shifts;
                if nd == 3
                    shifts_up = zeros([options.d1,options.d2,options.d3,3]);
                    do = size(shifts_up)./size(shifts_temp);
                    ds = size(shifts_temp);
                    dim = [options.d1,options.d2,options.d3];
                    if (0)
                        [Xq,Yq,Zq] = meshgrid(linspace(1,ds(2),dim(2)),linspace(1,ds(1),dim(1)),linspace(1,ds(3),dim(3)));
                    else
                        [Xq,Yq,Zq] = meshgrid(linspace((1+1/do(2))/2,ds(2)+(1-1/do(2))/2,dim(2)),linspace((1+1/do(1))/2,ds(1)+(1-1/do(1))/2,dim(1)),linspace((1+1/do(3))/2,ds(3)+(1-1/do(3))/2,dim(3)));
                        Xq(Xq<1)=1; Xq(Xq>dim(2))=dim(2);
                        Yq(Yq<1)=1; Yq(Yq>dim(1))=dim(1);
                        Zq(Zq<1)=1; Zq(Zq>dim(3))=dim(3);
                    end
                    if numel(shifts_temp) > 3
                        for dm = 1:3; shifts_up(:,:,:,dm) = interp3(shifts_temp(:,:,:,dm),Xq,Yq,Zq,'makima'); end
                        %tform = affine3d(diag([do([2,1,3])';1]));
                        %for dm = 1:3; shifts_up(:,:,:,dm) = imwarp(shifts_temp(:,:,:,dm),tform,'OutputView',imref3d([options.d1,options.d2,options.d3]),'SmoothEdges',true); end
                    else
                        for dm = 1:3; shifts_up(:,:,:,dm) = shifts_temp(dm); end
                    end
                    shifts_up(2:2:end,:,:,2) = shifts_up(2:2:end,:,:,2) + col_shift;
                    Mf = imwarp(Yt,-cat(4,shifts_up(:,:,:,2),shifts_up(:,:,:,1),shifts_up(:,:,:,3)),options.shifts_method,'FillValues',fill_value);
                else
                    shifts_up = imresize(shifts_temp,[options.d1,options.d2]);
                    shifts_up(2:2:end,:,2) = shifts_up(2:2:end,:,2) + col_shift;
                    Mf = imwarp(Yt,-cat(3,shifts_up(:,:,2),shifts_up(:,:,1)),options.shifts_method,'FillValues',fill_value);
                end

             Mf(Mf<minY) = minY;
             Mf(Mf>maxY) = maxY;
        end

        if ~strcmpi(options.output_type,'mat')
            rem_mem = rem(t,options.mem_batch_size);
            if rem_mem == 0; rem_mem = options.mem_batch_size; end
            if nd == 2; mem_buffer(:,:,rem_mem) = cast(Mf,data_type); end
            if nd == 3; mem_buffer(:,:,:,rem_mem) = cast(Mf,data_type); end
        end
        switch lower(options.output_type)
            case 'mat'
                if nd == 2; M_final(:,:,t) = cast(Mf,data_type); end
                if nd == 3; M_final(:,:,:,t) = cast(Mf,data_type); end
            case 'memmap'
                if rem_mem == options.mem_batch_size || t == T
                    if nd == 2; M_final.Y(:,:,t-rem_mem+1:t) = mem_buffer(:,:,1:rem_mem); end
                    if nd == 3; M_final.Y(:,:,:,t-rem_mem+1:t) = mem_buffer(:,:,:,1:rem_mem); end
                    M_final.Yr(:,t-rem_mem+1:t) = reshape(mem_buffer(1:d1*d2*d3*rem_mem),d1*d2*d3,rem_mem);
                end
            case {'hdf5','h5'}
                if rem_mem == options.mem_batch_size || t == T
                    if nd == 2; h5write(options.h5_filename,['/',options.h5_groupname],mem_buffer(:,:,1:rem_mem),[ones(1,nd),t-rem_mem+1],[sizY(1:nd),rem_mem]); end
                    if nd == 3; h5write(options.h5_filename,['/',options.h5_groupname],mem_buffer(:,:,:,1:rem_mem),[ones(1,nd),t-rem_mem+1],[sizY(1:nd),rem_mem]); end
                end
            case {'tif','tiff'}
                if rem_mem == options.mem_batch_size || t == T
                    saveastiff(cast(mem_buffer(:,:,1:rem_mem),data_type),options.tiff_filename,opts_tiff);
                end
        end

        if mod(t,bin_width) == 0 && upd_template
            if print_msg
                str=[num2str(t), ' out of ', num2str(T), ' frames registered, iteration ', num2str(it), ' out of ', num2str(iter), '..'];
                refreshdisp(str, prevstr, t);
                prevstr=str;
                %fprintf('%i out of %i frames registered, iteration %i out of %i \n',t,T,it,iter)
            end
            cnt_buf = cnt_buf + 1;
            if strcmpi(method{2},'mean')
                new_temp = cellfun(@(x) nanmean(x,nd+1), buffer, 'UniformOutput',false);
            elseif strcmpi(method{2},'median');
                new_temp = cellfun(@(x) nanmedian(x,nd+1), buffer, 'UniformOutput', false);
            end
            if strcmpi(method{1},'mean')
                cnt = t/bin_width + 1;
                template = cellfun(@plus, cellfun(@(x) x*(cnt-1)/cnt, template,'un',0), cellfun(@(x) x*1/cnt, new_temp,'un',0), 'un',0);
            elseif strcmpi(method{1},'median');
                if cnt_buf <= buffer_width
                    if nd == 2; buffer_med(:,:,cnt_buf) = cell2mat_ov(new_temp,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end
                    if nd == 3; buffer_med(:,:,:,cnt_buf) = cell2mat_ov(new_temp,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end
                else
                    buffer_med = circshift(buffer_med,[zeros(1,nd),-1]);
                    if nd == 2; buffer_med(:,:,buffer_width) = cell2mat_ov(new_temp,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end
                    if nd == 3; buffer_med(:,:,:,buffer_width) = cell2mat_ov(new_temp,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY); end
                end
                template = mat2cell_ov(nanmedian(buffer_med,nd+1),xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY);
            end
            fftTemp = cellfun(@fftn, template, 'un',0);
            temp_mat = cell2mat_ov(template,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY);
            fftTempMat = fftn(temp_mat);
        end

        if plot_flag && mod(t,1) == 0
            subplot(221); imagesc(Yt-add_value,[nnY,mmY]); title('Raw data','fontweight','bold','fontsize',14);
                            xlabel(sprintf('Frame %i out of %i',t,T),'fontweight','bold','fontsize',14); set(gca,'Xtick',[],'Ytick',[]);
            subplot(222); imagesc(Mf,[nnY,mmY]);  title('Motion Corrected','fontweight','bold','fontsize',14); colormap('bone'); axis off;
            subplot(223); quiver(shifts_up(:,:,:,1),shifts_up(:,:,:,2),'Autoscale','off'); title('Motion vector field','fontweight','bold','fontsize',14); axis off;
            subplot(224); imagesc(cell2mat_ov(template,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY)-add_value,[nnY,mmY]); title('Matching Template','fontweight','bold','fontsize',14); axis off
            drawnow;
            if make_avi
                currFrame = getframe(fig);
                writeVideo(vidObj,currFrame);
            end
        end
    end

if print_msg; fprintf('\n'); end

if it == iter
    template = cellfun(@(x) x - add_value,template,'un',0);
    template = cell2mat_ov(template,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap_pre,sizY);
end
if memmap
    M_final.shifts = shifts;
    M_final.template = template;
end

if make_avi && plot_flag
    close(vidObj);
end
maxNumCompThreads('automatic');
if print_msg; fprintf('done. \n'); end
end

author: flatironinstitute

repo: NoRMCorre


text:
% GCD    Greatest common divisor.
% G = GCD(A,B) is the greatest common divisor of A and B
function g = gcd(a,b)
	if isscalar(a)
		a = a*ones(size(b),'int64');
	elseif isscalar(b)
		b = a*ones(size(b),'int64');
	end
	if ~isequal(size(a),size(b))
		error('Size mismatch');
	end

	g = zeros(size(a),'int64');
	for k = 1:length(a)
		g(k) = a(k);
		while b(k) ~= 0
		   t    = b(k);
		   b(k) = mod(g(k),b(k));
		   g(k) = t;
		end
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function [clusternum, label] = struG2la(Z)
    [n,m] = size(Z);
    SS0 = sparse(n+m,n+m); SS0(1:n,n+1:end) = Z; SS0(n+1:end,1:n) = Z';
    [clusternum, label] = graphconncomp(SS0);
    label = label(1:n)';

end

author: CHLWR


text:
function unicycle_print ( n, u, title )

%*****************************************************************************80
%
%% UNICYCLE_PRINT prints a unicycle given in sequence form.
%
%  Example:
%
%    Input:
%
%      U = 7 1 4 5 2 3 6
%
%    Printed output:
%
%      7 1 4 5 2 3 6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the unicycle.
%
%    Input, integer U(N), the unicycle sequence vector.
%
%    Input, string TITLE, a title.
%
  inc = 20;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
    fprintf ( 1, '\n' );
  end

  for ilo = 1 : inc : n
    ihi = min ( n, ilo + inc - 1 );
    for i = ilo : ihi
      fprintf ( 1, '  %4d', u(i) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function I_out = vesselFilter(I_in, grid_spacing, scales, varargin)
%VESSELFILTER   Frangi's 3D vessel filter.
%
% DESCRIPTION:
%       vesselFilter filters a 3D image using Frangi's vessel filtering
%       algorithm [1-3]. The algorithm works by calculating the Hessian
%       matrix (containing second order gradients) at each image voxel. The
%       eigenvalues of this matrix are then ordered and used to classify
%       whether the voxel is part of a vessel. The algorithm is performed
%       over multiple scales by convolving the input image with a Gaussian
%       of the given input scales. The final output is taken as the maximum
%       of the vessel filtered image across all scales.
%
% USAGE:
%       I_out = vesselFilter(I_in, grid_spacing, scales)
%       I_out = vesselFilter(I_in, grid_spacing, scales, ...)
%
% INPUTS:
%       I_in         - 3D image data to filter
%       grid_spacing - [dx, dy, dz] in [m]
%       scales       - array of scales to use [m]
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       alpha        - Value of sensitivity parameter for metric that
%                      distinguishes between plate-like and other
%                      structures (vessel-like or ball-like).
%                      (default = 0.5)
%       beta         - Value of sensitivity parameter for metric that
%                      distinguishes between ball-like and other structures
%                      (vessel-like or plate-like).
%                      (default = 0.5)
%       c            - Value used to scale the sensitivity parameter for
%                      the noise metric (the sensitivity parameter itself
%                      is calculated automatically based on the magnitudes
%                      of the eigenvalues).
%                      (default = 1)
%       gamma        - normalisation factor for scale-space derivatives
%                      (default = 1)
%       DisplayUpdates - Boolean controlling whether command line updates
%                      are displayed (default = true)
%       Plot         - Boolean controlling whether a maximum intensity
%                      projection of the vessel filtered image at each
%                      scale is displayed (default = false)
%       Colormap     - colormap to use if 'Plot' is set to true (default =
%                      flipud(gray))
%
% OUTPUTS:
%       I_out        - vessel filtered image
%
% ABOUT:
%       author       - Bradley Treeby and Tanmayi Oruganti
%       date         - 9th May 2012
%       last update  - 21st August 2014
%
%       This function is based on original code by Dean Barratt and Yipeng
%       Hu, CMIC, UCL, 2006-2009
%
% REFERENCES:
%   [1] A. F. Frangi, W. J. Niessen, K. L. Vincken, M. A. Viergever (1998)
%       "Multiscale vessel enhancement filtering," MICCAI, pp. 130-137.
%   [2] R. Manniesing, M.A. Viergever, W.J. Niessen (2006) "Vessel
%       enhancing diffusion: A scale space representation of vessel
%       structures," Med. Image Anal. 10, pp. 815-25.
%   [3] T. Oruganti, J. Laufer, B. E. Treeby (2013) "Vessel filtering of
%       photoacoustic images," Proc. of SPIE, vol. 8581, p. 85811W-1.
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also interpftn, smooth

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% start timer
start_time = clock;

% set default parameters
num_req_input_variables = 3;    % minimum number of input variables
alpha = 0.5;                    % value of sensitivity parameter for plate vs vessel/ball measure
beta = 0.5;                     % value of sensitivity parameter for ball vs vessel/plate measure
c_scale = 1;                    % scale value of sensitivity parameter for noise
gamma = 0;                      % normalisation factor for scale-space derivatives
origin_smoothness_const = 1e-6; % constant for Manniesing's origin smoothness term
disp_updates = true;            % boolean controlling display of command line updates
plot_updates = false;           % boolean controlling MIP display
color_map = flipud(gray);       % colormap to use if images are displayed
cardan_tol = 0.001;             % tolerance for finding complex roots using cardanRoots

% check user defined inputs
% replace with user defined values if provided
if nargin < num_req_input_variables
    error('Incorrect number of inputs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'alpha'
                alpha = varargin{input_index + 1};
            case 'beta'
                beta = varargin{input_index + 1};
            case 'c'
                c_scale = varargin{input_index + 1};
            case 'gamma'
                gamma = varargin{input_index + 1};
            case 'Display'
                disp_updates = varargin{input_index + 1};
            case 'DisplayUpdates'
                disp_updates = varargin{input_index + 1};
            case 'Plot'
                plot_updates = varargin{input_index + 1};
            case 'Colormap'
                color_map = varargin{input_index + 1};
            otherwise
                error('Unknown optional input');
        end
    end
end

% get size of input volume
sz = size(I_in);

% create k-space grid to use for spectral derivatives
kgrid = makeGrid(sz(1), grid_spacing(1), sz(2), grid_spacing(2), sz(3), grid_spacing(3));

% force wavenumber vectors to be in the correct directions for bsxfun
kx_vec = kgrid.kx_vec;
ky_vec = kgrid.ky_vec.';
kz_vec = permute(kgrid.kz_vec, [2 3 1]);

% compute 3D fft of input image
Ik = fftn(I_in);

% preallocate empty output image
I_out = zeros(size(I_in));

% loop through the scales defined by the user
for scale_index = 1:length(scales)

    % extract the current scale
    scale = scales(scale_index);

    % update command line status
    if disp_updates
        disp(['Computing Vesselness filter for scale ' num2str(scale) ':']);
        tic, fprintf('  Calculating derivatives... ');
    end

    % create normalised frequency domain gaussian at current scale
    FG = exp( (-0.5*scale^2)*(kgrid.kx.^2 + kgrid.ky.^2 + kgrid.kz.^2) );

    % compute components of Hessian matrix convolved with Gaussian
    Hxx = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*kx_vec).^2, FG)) )), [], 1);
    Hyy = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*ky_vec).^2, FG)) )), [], 1);
    Hzz = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*kz_vec).^2, FG)) )), [], 1);
    Hxy = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*kx_vec), bsxfun(@times, (1i*ky_vec), FG))) )), [], 1);
    Hxz = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*kx_vec), bsxfun(@times, (1i*kz_vec), FG))) )), [], 1);
    Hyz = reshape(real(ifftn( Ik .* ifftshift( bsxfun(@times, (1i*ky_vec), bsxfun(@times, (1i*kz_vec), FG))) )), [], 1);

    % update command line status
    if disp_updates
        toc, tic, fprintf('  Computing eigenvalues... ');
    end

    % compute trace of Hessian matrix
    P2 = -(Hxx + Hyy + Hzz);

    % compute principal minors of Hessian matrix
    M11 = Hyy.*Hzz - Hyz.^2;
    M22 = Hzz.*Hxx - Hxz.^2;
    M33 = Hxx.*Hyy - Hxy.^2;
    P1 = (M11 + M22 + M33);

    % compute determinant of Hessian matrix
    P0 = - Hxx.*M11 ...
         + Hxy.*(Hxy.*Hzz - Hyz.*Hxz) ...
         - Hxz.*(Hxy.*Hyz - Hyy.*Hxz);

    % find eigenvalues using roots of 3rd order polynomial
    P3 = 1;
    eigenval_mat = cardanRoots(P3, P2, P1, P0, cardan_tol).';

    % update command line status
    if disp_updates
        toc, tic, fprintf('  Sorting eigenvalues... ');
    end

    % sort eigenvalues based on their absolute values
    [~, ind] = sort(abs(eigenval_mat));
    ind = bsxfun(@plus, (0:numel(I_in)-1)*3, ind);
    eigenval_mat = eigenval_mat(ind);

    % update command line status
    if disp_updates
        toc, tic, fprintf('  Evaluating vesselness measure... ');
    end

    % calculate noise parameter dynamically if not given
    Rs = sqrt( abs(eigenval_mat(1, :)).^2 + abs(eigenval_mat(2, :)).^2 + abs(eigenval_mat(3, :)).^2 );
    c = c_scale*0.5*max(Rs(:));

    % calculate overall vesselness measure and scale by Lindeberg?s constant
    %   term 1: plate vs vessel/ball measure
    %   term 2: ball vs vessel/plate measure
    %   term 3: noise measure (Frobenius matrix norm of the Hessian)
    %   term 4: Manniesing's origin smoothness term
    V = scale.^(gamma) .* ...
        (1 - exp(-(abs(eigenval_mat(2, :))./abs(eigenval_mat(3, :))).^2./(2*alpha^2))) .* ...
        (exp(-(abs(eigenval_mat(1, :))./sqrt(abs(eigenval_mat(2, :).*eigenval_mat(3, :)))).^2./(2*beta^2))) .* ...
        (1 - exp(-Rs.^2./(2*c^2))) .* ...
        (exp(-(2*origin_smoothness_const.^2)./(abs(eigenval_mat(2, :)).*eigenval_mat(3, :).^2)));

    % eliminate structures with positive eigenvalues as these indicate
    % regions of low optical absorption on a background of high optical
    % absorption
    V(eigenval_mat(2, :) > 0 | eigenval_mat(3, :) > 0) = 0;

    % reshape to be the correct size
    V = reshape(V, sz);

    % take maximum response over multiple scales
    I_out = max(I_out, V);

    % display time elapsed
    if disp_updates
        toc
    end

    % plot steps if required
    if plot_updates

        % get suitable axis scaling
        [~, axis_scale, axis_prefix] = scaleSI(max([sz(1)*grid_spacing(1), sz(2)*grid_spacing(2)]));

        % plot figure
        figure;
        imagesc(kgrid.y_vec * axis_scale, kgrid.x_vec * axis_scale,  max(V, [], 3));
        axis image;
        colormap(color_map);
        title(['Vessel Filtered Image (MIP), Scale = ' num2str(scale)]);
        ylabel(['x [' axis_prefix 'm]']);
        xlabel(['y [' axis_prefix 'm]']);
        colorbar;
        drawnow;

        % plot final multi-scale image
        if scale_index == length(scales)
            figure;
            imagesc(kgrid.y_vec * axis_scale, kgrid.x_vec * axis_scale,  max(I_out, [], 3));
            axis image;
            colormap(color_map);
            title('Vessel Filtered Image (MIP), Multi-Scale');
            ylabel(['x [' axis_prefix 'm]']);
            xlabel(['y [' axis_prefix 'm]']);
            colorbar;
            drawnow;
        end
    end
end

% display total time elapsed
if disp_updates
    disp(['  Total computation time ' scaleTime(etime(clock, start_time))]);
end

% end of function
end

% =========================================================================
% SUB FUNCTIONS
% =========================================================================

function roots = cardanRoots(varargin)
% Find roots of third order polynomials using Cardan's formula
%
% INPUT
%   P: (n x 4) array, each row corresponds to coefficients of each
%   polynomial, P(:,1)*x^3 + P(:,2)*x^2 + P(:,3)*x + P(:,4)
% OUTPUT
%   roots: (n x 3) array, each row correspond to the roots of P
%
% To adjust the parameter below which the the discriminant is considerered
% as nil, use
%   CardanRoots(P, tol)
% Adjusting tol is useful to avoid the real roots become complex due to
% numerical accuracy. The default TOL is 0
%
% http://www.sosmath.com/algebra/factor/fac11/fac11.html
% http://mathforum.org/dr.math/faq/faq.cubic.equations.html
%
% See also: roots, ParabolaRoots, eig3
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History:
%     Original 20-May-2010

% Copyright (c) 2010, Bruno Luong
% 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
%
% 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.

% Adjustable parameter
tol = 0;

if nargin<4
    P = varargin{1};
    a = P(:,1);
    b = P(:,2);
    c = P(:,3);
    d = P(:,4);
    if nargin>=2
        tol = varargin{2};
    end
else
    [a, b, c, d] = deal(varargin{1:4});
    if nargin>=5
        tol = varargin{5};
    end
end

if ~isequal(a,1)
    b = b./a;
    c = c./a;
    d = d./a;
end

b2 = b.^2;

p = -b2/3 + c;
q = ((2/27)*b2-(1/3)*c).*b + d;
delta = q.^2 + (4/27)*p.^3;

% Three cases of discriminant sign
iscmplx = imag(p) | imag(q);
notcmplx = ~iscmplx;
deltanull = notcmplx & abs(delta)<tol; % = 0
deltaneg = notcmplx & delta<0;
deltapos = notcmplx & ~(deltanull | deltaneg);

n = size(delta,1);
roots = zeros(n, 3, class(delta));

if any(deltanull)
    idx = find(deltanull);
    roots(idx,:) = CardanNull(p(idx), q(idx));
end

if any(deltaneg)
    idx = find(deltaneg);
    roots(idx,:) = CardanNeg(q(idx), delta(idx));
end

if any(deltapos)
    idx = find(deltapos);
    roots(idx,:) = CardanPos(q(idx), delta(idx));
end

if any(iscmplx)
    idx = find(iscmplx);
    roots(idx,:) = CardanCmplx(p(idx), q(idx), delta(idx));
end

roots = bsxfun(@minus, roots, b/3);

end

function roots = CardanNull(p, q)

    S1 = 3*q./p;
    S2 = -0.5*S1;
    roots = [S1 S2 S2];  % double real solutions

end

function roots = CardanNeg(q, delta)

    alfa = -q;
    beta = sqrt(-delta);
    r2 = alfa.^2-delta;
    rho = (4^(1/3))*exp(log(r2)/6);
    theta = atan2(beta,alfa)/3;
    alfa = rho.*cos(theta);
    beta = rho.*sin(theta);
    S1 = alfa;
    x = (-0.5)*alfa;
    y = (sqrt(3)/2)*beta;
    S2 = x-y;
    S3 = x+y;
    roots = [S1 S2 S3];

end

function roots = CardanPos(q, delta)

    sqrtdelta = sqrt(delta);
    u3 = (-q+sqrtdelta)/2;
    v3 = (-q-sqrtdelta)/2;

    % Cubic roots of u3 and v3
    u = sign(u3).*exp(log(abs(u3))/3);
    v = sign(v3).*exp(log(abs(v3))/3);

    S1 = u+v;
    % Complex solutions
    j = complex(-0.5,sqrt(3)/2);
    j2 = complex(-0.5,-sqrt(3)/2);
    S2 = j*u+j2*v;
    S3 = conj(S2);
    roots = [S1 S2 S3];

end

function roots = CardanCmplx(p, q, delta)

    sqrtdelta = sqrt(delta);
    u3 = (-q+sqrtdelta)/2;
    v3 = (-q-sqrtdelta)/2;

    % we need u*v = -p/3
    p = (-1/3)*p;
    iu = abs(u3)>abs(v3);
    u = zeros(size(u3),class(u3));
    v = zeros(size(v3),class(v3));
    u(iu) = exp(log(u3(iu))/3);
    v(iu) = p(iu)./u(iu);
    v(~iu) = exp(log(v3(~iu))/3);
    u(~iu) = p(~iu)./v(~iu);

    S1 = u+v;

    j = complex(-0.5,sqrt(3)/2);
    j2 = complex(-0.5,-sqrt(3)/2);
    S2 = j*u+j2*v;
    S3 = j2*u+j*v;
    roots = [S1 S2 S3];

end

repo: Aero-matlab


text:
function Decis = ErrEval (FuncName, ErrCode);
%
%	Decis = ErrEval (FuncName, ErrCode);
%
%This function displays an error message based on ErrCode and
%returns 1 if it's a lethal error (Err_) , or 0 if it's
%an error or a warning (Wrn_)that does not require aborting a function
%
%
%ErrCode is a string made out of two parts
%  The first part is a 4 letter string containing
%   Err_ or Wrn_  depending on wether it's an error or
%   a warning.
%  The second part is a string that indicates the type of
%   error encountered.
%
%Defined values of ErrCode are :
%  No						%No problem
%  BadInpSize			%One of the input variables has a bad size
%  BadInOut          %Bad combination of input and output parameters
%  FewInp				%Not enough inputs
%  ManyInp				%Too many inputs
%  FewOut				%Not enough outputs
%  ManyOut				%Too many outputs
%  SubFunc				%An error has occurred in a function within the script
%  BadOpt            %Bad Option
%  BadOptCombo       %Bad Option Combination
%  InpFieldMiss      %A field is missing from one of the input structures
%  OptFieldMiss      %A fiels is missing from the Options structures
%  EmptyInp          %Input data is empty
%  FileExist         %File exists
%  FileNotExist      %File does not exist
% To pass a generic message, replace the ErrCode by any ther string
%
% Here's a typical example of the function usage
% if (ErrEval ('LoadSliceData','Err_Image Type Not Supported')), err = 1;	return;end
%
%
%		Ziad Saad Thu Mar 26 11:43:51 CST 1998

if (eq_str(ErrCode,'No')), %in case you send a No only to
	Decis = 0;
	return;
end

if (length(ErrCode) < 5),
	fprintf (2,'\nError in ErrEval, ErrCode %s is ambiguous.\nReturning a 1 decision.\n\n\a',ErrCode);
	Decis = 1;
	return;
end

tmp = ErrCode(1:3);

switch tmp
	case 'Err'
		Decis = 1;
	case 'Wrn'
		Decis = 0;
	otherwise
		warning ('Error in ErrEval, prefix %s from %s is ambiguous.\nReturning a 1 decision.\n\n\a',tmp, ErrCode);
		Decis = 1;
		return;
end

tmp = ErrCode(5:length(ErrCode));

switch tmp
	case 'No'
		Decis = 0;		%even if it was Err_No
		return;
	case 'BadInpSize'
		s = 'Bad size of at least one of the input parameters.';
	case 'BadInOut',
		s = 'Bad combination of input and output parameters.';
	case 'FewInp'
		s = 'Too few inputs.';
	case 'ManyInp'
		s = 'Too many inputs.';
	case 'FewOut'
		s = 'Too few outputs';
	case 'ManyOut'
		s = 'Too many outputs';
	case 'SubFunc'
		s = 'Error (or Warning) in a function within this function (or script).';
	case 'BadOpt'
		s = 'Bad Option';
	case 'InpFieldMiss',
		s = 'A field is missing from one of the input structures.';
	case 'OptFieldMiss',
		s = 'A field is missing from the Options structure.';
	case 'BadOptCombo',
		s = 'Bad Combinations of Options.';
	case 'EmptyInp',
		s = 'One or all of input data is empty.';
	case 'FileExist',
		s = 'File specified exists.';
	case 'FileNotExist',
		s = 'File specified does not exist.';
	otherwise
		s = tmp;
end

%For the error message
if (Decis),
	warning('\a\nError in %s : %s\n\n',FuncName,s);
else
	warning('\nWarning from %s : %s\n\n',FuncName,s);
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function lattice_rule_test085 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST085 tests LATTICE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994, page 18.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST085\n' );
  fprintf ( 1, '  LATTICE is a lattice rule for periodic functions.\n' );
  fprintf ( 1, '  However, we apply it to a nonperiodic function\n' );
  fprintf ( 1, '  just to see how it does.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );

  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  exact = e_01_2d ( dim_num, a, b );

  z(1:dim_num) = [ 1, 2 ];
  i4vec_print ( dim_num, z, '  The lattice generator vector:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         K         M    EXACT       ESTIMATE    ERROR\n' );
  fprintf ( 1, '\n' );

  for k = 3 : 18

    m = fibonacci ( k );

    quad = lattice ( dim_num, m, z, @f_01_2d );

    error = abs ( exact - quad );

    fprintf ( 1, '  %8d  %8d  %10.6f  %10.6f  %10.6e\n', k, m, exact, quad, error );

  end

  return
end


text:
% readSPE.m
function image = readSPE(dirPath, filename)
% readSPE.m  Read princeton instrument *.SPE image file
%
% image = readSPE(filename)
%       = readSPE(dirPath, filename)
%
% where: image - 3D array of image(s)
%     filename - path and filename (string)
%      dirPath - optional directory path (string)
%
% Image is returned as a 3D array, where each image is the first 2
% dimensions, and successive images are stored along the 3rd dimension.
%
% The image is stored as it is shown in WinVIEW; the first two dimensions
% are stored as [pixel,stripe]
%
% returns -1 if unsuccessful, though if any of the low level access
% functions error out, this method will just fail with their error
% descriptions. Only the uint16 pixel datatype feature has been tested, so
% be more cautious than usual if using other datatypes.
%
% There is far more information stored in the SPE header (4100 bytes).
% This simple method only pulls the information that is required for the
% author's data processing. For further documentation of the header byte
% structure, seek the sources I used:
%
%  Matt Clarke's Python Script Description:
%   http://www.scipy.org/Cookbook/Reading_SPE_files
%
%  Stuart B. Wilkins's Python Code:
%   https://pyspec.svn.sourceforge.net/svnroot/pyspec/trunk/pyspec/ccd/files.py
%
%  ImageJ's plugin to import SPE files, written in Java:
%   http://rsbweb.nih.gov/ij/
%

% Author: Carl Hall
% Date: June 2012

% Modification History
%  April 2012 - Original Code
%   June 2012 - Modified to maintain output array as same datatype as SPE
%               file

%% SPE Header byte Structure (not fully checked)
%
% Offset    Bytes    Type     Description
% 0x0000    2        int      Controller Version
% 0x0002    2        int      Logic Output
% 0x0004    2        int      AppHiCapLowNoise
% 0x0006    2        int      dxDim
% 0x0008    2        int      timingMode
% 0x000A    4        float    Exposure
% 0x000E    2        int      vxDim
% 0x0010    2        int      vyDim
% 0x0012    2        int      dyDim
% 0x0014    10       string   Date
% 0x0024    4        float    DetTemperature
% 0x0028    2        int      DectectorType
% 0x002A    2        int      xDim
% 0x002C    2        int      TriggerDiode
% 0x002E    4        float    DelayTime
% 0x0032    2        int      ShutterControl
% 0x0034    2        int      AbsorbLive
% 0x0036    2        int      AbsorbMode
% 0x0038    2        int      canDoVirtualChip
% 0x003A    2        int      thresholdMinLive
% 0x003C    4        float    threshholdMin
% 0x0040    2        int      thresholdMaxLive
% 0x0042    4        float    threshholdMax
% 0x006C    2        int      Data Type
% 0x00AC    7        string   Time
% 0x00BC    2        int      ADCOffset
% 0x00BE    2        int      ADCRate
% 0x00C0    2        int      ADCType
% 0x00C2    2        int      ADCRes
% 0x00C4    2        int      ADCBitAdj
% 0x00C6    2        int      Gain
% 0x00C8    80       string   Comments
% 0x0118    80       string   Comments
% 0x0168    80       string   Comments
% 0x01B8    80       string   Comments
% 0x0208    80       string   Comments
% 0x0258    2        int      GeometricOps
% 0x0290    2        int      ydim
% 0x05A6    2        uint32   zdim
% 0x05D0    2        int      NumROIExperiment
% 0x05D2    2        int      NumROI
% 0x05D4    60       int[]    allROI
%
% 0x1004    ...      datatype   Image Data

%% Start of Code

% parse optional input
if nargin>1
    filename = strcat(dirPath,filename);
else
    filename = dirPath;
end

% Open the file
fd = fopen(filename,'r');
if(fd < 0)
    error('Could not open file, bad filename')
end

% Get the image dimensions:
stripDim = getData(fd, '2A', 'uint16');     %first dim
pixelDim = getData(fd, '290', 'uint16');    %second dim
nDim = getData(fd, '5A6', 'uint32');        %third dim

% Get the pixel data type
dataType = getData(fd, '6C', 'uint16');

% Get the image
fseek(fd, hex2dec('1004'), 'bof');

image = zeros([pixelDim,stripDim,nDim]);
switch dataType
    case 0     % single precision float (4 bytes)
        image = single(image);      %maintain datatype in function output
        for i=1:nDim
            image(:,:,i) = fread(fd, [stripDim,pixelDim], 'float32')';
        end

    case 1    % long int (4 bytes)
        image = int32(image);
        for i=1:nDim
            image(:,:,i) = fread(fd, [stripDim,pixelDim], 'int32')';
        end

    case 2    % short int (2 bytes)
        image = int16(image);
        for i=1:nDim
            image(:,:,i) = fread(fd, [stripDim,pixelDim], 'int16')';
        end

    case 3    % short unsigned int (2 bytes)
        image = uint16(image);
        for i=1:nDim
            image(:,:,i) = fread(fd, [stripDim,pixelDim], 'uint16')';
        end

    otherwise
        image = -1;
end
end

%%
% getData() reads one piece of data at a specific location
%
function data = getData(fd, hexLoc, dataType)
% Inputs: fd - int    - file descriptor
%     hexLoc - string - location of data relative to beginning of file
%   dataType - string - type of data to be read
%
fseek(fd, hex2dec(hexLoc), 'bof');
data = fread(fd, 1, dataType);
end


text:
%TRPLOT Draw a coordinate frame
%
% TRPLOT(T, OPTIONS) draws a 3D coordinate frame represented by the homogeneous
% transform T (4x4).
%
% H = TRPLOT(T, OPTIONS) as above but returns a handle.
%
% TRPLOT(H, T) moves the coordinate frame described by the handle H to
% the pose T (4x4).
%
% TRPLOT(R, OPTIONS) as above but the coordinate frame is rotated about the
% origin according to the orthonormal rotation matrix R (3x3).
%
% H = TRPLOT(R, OPTIONS) as above but returns a handle.
%
% TRPLOT(H, R) moves the coordinate frame described by the handle H to
% the orientation R.
%
% Options::
% 'color',C          The color to draw the axes, MATLAB colorspec C
% 'noaxes'           Don't display axes on the plot
% 'axis',A           Set dimensions of the MATLAB axes to A=[xmin xmax ymin ymax zmin zmax]
% 'frame',F          The coordinate frame is named {F} and the subscript on the axis labels is F.
% 'text_opts', opt   A cell array of MATLAB text properties
% 'handle',H         Draw in the MATLAB axes specified by the axis handle H
% 'view',V           Set plot view parameters V=[az el] angles, or 'auto'
%                    for view toward origin of coordinate frame
% 'length',s         Length of the coordinate frame arms (default 1)
% 'arrow'            Use arrows rather than line segments for the axes
% 'width', w         Width of arrow tips (default 1)
% 'thick',t          Thickness of lines (default 0.5)
% '3d'               Plot in 3D using anaglyph graphics
% 'anaglyph',A       Specify anaglyph colors for '3d' as 2 characters for
%                    left and right (default colors 'rc'): chosen from
%                    r)ed, g)reen, b)lue, c)yan, m)agenta.
% 'dispar',D         Disparity for 3d display (default 0.1)
% 'text'             Enable display of X,Y,Z labels on the frame
% 'labels',L         Label the X,Y,Z axes with the 1st, 2nd, 3rd character of the string L
% 'rgb'              Display X,Y,Z axes in colors red, green, blue respectively
% 'rviz'             Display chunky rviz style axes
%
% Examples::
%
%       trplot(T, 'frame', 'A')
%       trplot(T, 'frame', 'A', 'color', 'b')
%       trplot(T1, 'frame', 'A', 'text_opts', {'FontSize', 10, 'FontWeight', 'bold'})
%       trplot(T1, 'labels', 'NOA');
%
%       h = trplot(T, 'frame', 'A', 'color', 'b');
%       trplot(h, T2);
%
% 3D anaglyph plot
%       trplot(T, '3d');
%
% Notes::
% - The 'rviz' option is equivalent to 'rgb', 'notext', 'noarrow',
%   'thick', 5.
% - The arrow option requires the third party package arrow3 from File
%   Exchange.
% - The handle H is an hgtransform object.
% - When using the form TRPLOT(H, ...) to animate a frame it is best to set
%   the axis bounds.
% - The '3d' option requires that the plot is viewed with anaglyph glasses.
% - You cannot specify 'color' and '3d' at the same time.
%
% See also TRPLOT2, TRANIMATE.

%TODO:
% 'rviz', chunky RGB lines, no arrows

% Copyright (C) 1993-2015, 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

% TODO
%  need to decide how to handle scaling
%  what does hold on mean?  don't touch scaling?

function hout = trplot(T, varargin)

    if isscalar(T) && ishandle(T)
        % trplot(H, T)
        H = T; T = varargin{1};
        if isrot(T)
            T = r2t(T);
        end
        set(H, 'Matrix', T);

        % for the 3D case retrieve the right hgtransform and set it
        hg2 = get(H, 'UserData');
        if ~isempty(hg2)
            set(hg2, 'Matrix', T);
        end

        return;
    end

    if size(T,3) > 1
        error('trplot cannot operate on a sequence');
    end
    if ~ishomog(T) && ~isrot(T)
        error('trplot operates only on transform (4x4) or rotation matrix (3x3)');
    end

    opt.color = [];
    opt.rgb = false;
    opt.axes = true;
    opt.axis = [];
    opt.frame = [];
    opt.text_opts = [];
    opt.view = [];
    opt.width = 1;
    opt.arrow = false;
    opt.labels = 'XYZ';
    opt.handle = [];
    opt.anaglyph = 'rc';
    opt.d_3d = false;
    opt.dispar = 0.1;
    opt.thick = 0.5;
    opt.length = 1;
    opt.text = true;
    opt.lefty = false;
    opt.rviz = false;

    opt = tb_optparse(opt, varargin);

    if opt.rviz
        opt.thick = 5;
        opt.arrow = false;
        opt.rgb = true;
        opt.text = false;
    end

    if ~isempty(opt.color) && opt.d_3d
        error('cannot specify ''color'' and ''3d'', use ''anaglyph'' option');
    end
    if isempty(opt.color)
        opt.color = 'b';
    end
    if isempty(opt.text_opts)
        opt.text_opts = {};
    end

    if opt.d_3d
        opt.color = ag_color(opt.anaglyph(1));
    end

    if isempty(opt.axis)
        % determine some default axis dimensions

        % get the origin of the frame
        if isrot(T)
            c = [0 0 0];  % at zero for a rotation matrix
        else
            c = transl(T);
        end

        d = 1.2;
        opt.axis = [c(1)-d c(1)+d c(2)-d c(2)+d c(3)-d c(3)+d];

    end

    % TODO: should do the 2D case as well

    if ~isempty(opt.handle)
        hax = opt.handle;
        hold(hax);
    else
        ih = ishold;
        if ~ih
            % if hold is not on, then clear the axes and set scaling
            cla
            if ~isempty(opt.axis)
                axis(opt.axis);
            end
            daspect([1 1 1]);

            if opt.axes
                xlabel( 'X');
                ylabel( 'Y');
                zlabel( 'Z');
                rotate3d on
            end
            new_plot = true;
        end
        hax = gca;
        hold on
    end
    % hax is the handle for the axis we will work with, either new or
    % passed by option 'handle'

    opt.text_opts = [opt.text_opts, 'Color', opt.color];

    hg = hgtransform('Parent', hax);

    % trplot( Q.R, fmt, color);
    if isrot(T)
        T = r2t(T);
    end

    % create unit vectors
    o =  [0 0 0]';
    x1 = opt.length*[1 0 0]';
    y1 = opt.length*[0 1 0]';
    if opt.lefty
        z1 = opt.length*[0 0 -1]';
    else
        z1 = opt.length*[0 0 1]';
    end

    % draw the axes

    mstart = [o o o]';
    mend = [x1 y1 z1]';

    if opt.rgb
        axcolors = {'r', 'g', 'b'};
    else
        axcolors = { opt.color, opt.color, opt.color};
    end

    if opt.arrow
%         % draw the 3 arrows
%         S = [opt.color num2str(opt.width)];
%         ha = arrow3(mstart, mend, S);
%         for h=ha'
%             set(h, 'Parent', hg);
%         end
          daspect([1,1,1])
          for i=1:3
              ha = arrow3(mstart(i,1:3), mend(i,1:3), [axcolors{i} num2str(opt.width)]);
              set(ha, 'Parent', hg);
          end
    else
        for i=1:3
            plot2([mstart(i,1:3); mend(i,1:3)], 'Color', axcolors{i}, ...
                'LineWidth', opt.thick, ...
                'Parent', hg);
        end
    end

    % label the axes
    if isempty(opt.frame)
        fmt = '%c';
    else
        fmt = sprintf('%%c_{%s}', opt.frame);
    end

    if opt.text
        % add the labels to each axis
        h = text(x1(1), x1(2), x1(3), sprintf(fmt, opt.labels(1)), 'Parent', hg);
        set(h, opt.text_opts{:});

        h = text(y1(1), y1(2), y1(3), sprintf(fmt, opt.labels(2)), 'Parent', hg);
        set(h, opt.text_opts{:});

        h = text(z1(1), z1(2), z1(3), sprintf(fmt, opt.labels(3)), 'Parent', hg);
        set(h, opt.text_opts{:});
    end

    % label the frame
    if ~isempty(opt.frame)
        h = text(o(1)-0.04*x1(1), o(2)-0.04*y1(2), o(3)-0.04*z1(3), ...
            ['\{' opt.frame '\}'], 'Parent', hg);
        set(h, 'VerticalAlignment', 'middle', ...
            'HorizontalAlignment', 'center', opt.text_opts{:});
    end

    if ~opt.axes
        set(gca, 'visible', 'off');
    end
    if ischar(opt.view) && strcmp(opt.view, 'auto')
        cam = x1+y1+z1;
        view(cam(1:3));
    elseif ~isempty(opt.view)
        view(opt.view);
    end
    if isempty(opt.handle) && ~ih
        grid on
        hold off
    end

    % now place the frame in the desired pose
    set(hg, 'Matrix', T);

    if opt.d_3d
        % 3D display.  The original axes are for the left eye, and we add
        % another set of axes to the figure for the right eye view and
        % displace its camera to the right of that of that for the left eye.
        % Then we recursively call trplot() to create the right eye view.

        left = gca;
        right = axes;

        % compute the offset in world coordinates
        off = -t2r(view(left))'*[opt.dispar 0 0]';
        pos = get(left, 'CameraPosition');

        set(right, 'CameraPosition', pos+off');
        set(right, 'CameraViewAngle', get(left, 'CameraViewAngle'));
        set(right, 'CameraUpVector', get(left, 'CameraUpVector'));
        target = get(left, 'CameraTarget');
        set(right, 'CameraTarget', target+off');

        % set perspective projections
                set(left, 'Projection', 'perspective');
        set(right, 'Projection', 'perspective');

        % turn off axes for right view
        set(right, 'Visible', 'Off');

        % set color for right view
        hg2 = trplot(T, 'color', ag_color(opt.anaglyph(2)));

        % the hgtransform for the right view is user data for the left
        % view hgtransform, we need to change both when we rotate the
        % frame.
        set(hg, 'UserData', hg2);
    end

    % optionally return the handle, for later modification of pose
    if nargout > 0
        hout = hg;
    end
end

function out = ag_color(c)

% map color character to an color triple, same as anaglyph.m

    % map single letter color codes to image planes
    switch c
    case 'r'
        out = [1 0 0];        % red
    case 'g'
        out = [0 1 0];        % green
    case 'b'
        % blue
        out = [0 0 1];
    case 'c'
        out = [0 1 1];        % cyan
    case 'm'
        out = [1 0 1];        % magenta
    case 'o'
        out = [1 1 0];        % orange
    end
end

author: Allopart

repo: rbpf-gmapping