function [file_name_to_save, filename_raw, filename_reg] = min1pipe(Fsi, Fsi_new, spatialr, se, ismc, flag) % main_processing % need to decide whether to use parallel computing % Fsi: raw sampling rate % Fsi_new: in use sampling rate % spatialr: spatial downsampling factor % Jinghao Lu 06/10/2016 %% configure paths %% min1pipe_init; %% initialize parameters %% defpar = default_parameters; aflag = false; if nargin < 1 || isempty(Fsi) Fsi = defpar.Fsi; end if nargin < 2 || isempty(Fsi_new) Fsi_new = defpar.Fsi_new; end if nargin < 3 || isempty(spatialr) spatialr = defpar.spatialr; aflag = true; end if nargin < 4 || isempty(se) se = defpar.neuron_size; aflag = true; end if nargin < 5 || isempty(ismc) ismc = true; end if nargin < 6 || isempty(flag) flag = 1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%% parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% user defined parameters %%% %%% Params.Fsi = Fsi; %%% Params.Fsi_new = Fsi_new; %%% Params.spatialr = spatialr; %%% Params.neuron_size = se; %%% half neuron size; 9 for Inscopix and 5 %%% %%% for UCLA, with 0.5 spatialr separately %%% %%% %%% fixed parameters (change not recommanded) %%% %%% Params.anidenoise_iter = 4; %%% denoise iteration %%% Params.anidenoise_dt = 1/7; %%% denoise step size %%% Params.anidenoise_kappa = 0.5; %%% denoise gradient threshold %%% Params.anidenoise_opt = 1; %%% denoise kernel choice %%% Params.anidenoise_ispara = 1; %%% if parallel (denoise) %%% Params.bg_remove_ispara = 1; %%% if parallel (backgrond removal) %%% Params.mc_scl = 0.004; %%% movement correction threshold scale %%% Params.mc_sigma_x = 5; %%% movement correction spatial uncertainty %%% Params.mc_sigma_f = 10; %%% movement correction fluid reg weight %%% Params.mc_sigma_d = 1; %%% movement correction diffusion reg weight %%% Params.pix_select_sigthres = 0.8; %%% seeds select signal level %%% Params.pix_select_corrthres = 0.6; %%% merge correlation threshold1 %%% Params.refine_roi_ispara = 1; %%% if parallel (refine roi) %%% Params.merge_roi_corrthres = 0.9; %%% merge correlation threshold2 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% get dataset info %% [path_name, file_base, file_fmt] = data_info; hpipe = tic; for i = 1: length(file_base) %%% judge whether do the processing %%% filecur = [path_name, file_base{i}, '_data_processed.mat']; msg = 'Redo the analysis? (y/n)'; overwrite_flag = judge_file(filecur, msg); if overwrite_flag %% data cat %% %%% --------- 1st section ---------- %%% Fsi = Params.Fsi; Fsi_new = Params.Fsi_new; spatialr = 1; [m, filename_raw, imaxn, imeanf, pixh, pixw, nf, imx1, imn1] = data_cat(path_name, file_base{i}, file_fmt{i}, Fsi, Fsi_new, spatialr); %%% remove dead pixels %%% [m, imaxn] = remove_dp(m, 'frame_allt'); %%% spatial downsampling after auto-detection %%% [m, Params, pixh, pixw] = downsamp(path_name, file_base{i}, m, Params, aflag, imaxn); %% neural enhancing batch version %% %%% --------- 2nd section ---------- %%% filename_reg = [path_name, file_base{i}, '_reg.mat']; [m, imaxy1, overwrite_flag, imx2, imn2, ibmean] = neural_enhance(m, filename_reg, Params); %% neural enhancing postprocess %% if overwrite_flag nflag = 1; m = noise_suppress(m, imaxy1, Fsi_new, nflag); end %% movement correction %% if ismc if overwrite_flag pixs = min(pixh, pixw); Params.mc_pixs = pixs; Fsi_new = Params.Fsi_new; scl = Params.neuron_size / (7 * pixs); sigma_x = Params.mc_sigma_x; sigma_f = Params.mc_sigma_f; sigma_d = Params.mc_sigma_d; se = Params.neuron_size; [m, corr_score, raw_score, scl, imaxy] = frame_reg(m, imaxy1, se, Fsi_new, pixs, scl, sigma_x, sigma_f, sigma_d); Params.mc_scl = scl; %%% update latest scl %%% % file_name_to_save = [path_name, file_base{i}, '_data_processed.mat']; % if exist(file_name_to_save, 'file') % delete(file_name_to_save) % end save(m.Properties.Source, 'corr_score', 'raw_score', 'imaxy', '-v7.3', '-append'); else imaxy = m.imaxy; end else if overwrite_flag m = frame_stab(m); %%% spatiotemporal stabilization %%% end imaxy = imaxy1; end %% movement correction postprocess %% %%% --------- 3rd section ---------- %%% nflag = 2; filename_reg_post = [path_name, file_base{i}, '_reg_post.mat']; m = noise_suppress(m, imaxy, Fsi_new, nflag, filename_reg_post); %% get rough roi domain %% mask = dominant_patch(imaxy); %% parameter init %% [P, options] = par_init(m); %% select pixel %% [sigrf, roirf, seedsupdt, bgrf, bgfrf, datasmthf1, cutofff1, pkcutofff1] = iter_seeds_select(m, mask, Params, P, options, flag); %% merge roi %% corrthres = Params.merge_roi_corrthres; [roimrg, sigmrg, seedsmrg, datasmthf2, cutofff2, pkcutofff2] = merge_roi(m, roirf, sigrf, seedsupdt, imaxy, datasmthf1, cutofff1, pkcutofff1, corrthres); % %% 2nd step clean seeds %% % sz = Params.neuron_size; % [roic, sigc, seedsc, datasmthc, cutoffc, pkcutoffc] = final_seeds_select(m, roimrg, sigmrg, seedsmrg, datasmthf2, cutofff2, pkcutofff2, sz, imax); %% refine roi again %% noise = P.sn; Puse.p = 0; Puse.options = options; Puse.noise = noise; ispara = Params.refine_roi_ispara; [roifn1, sigfn1, seedsfn1, datasmthfn1, cutofffn1, pkcutofffn1] = refine_roi(m, sigmrg, bgfrf, roimrg, seedsmrg, Puse.noise, datasmthf2, cutofff2, pkcutofff2, ispara); [bgfn, bgffn] = bg_update(m, roifn1, sigfn1); %% refine sig again %% Puse.p = 2; %%% 2nd ar model used %%% Puse.options.p = 2; Puse.options.temporal_iter = 1; [sigfn1, bgffn, roifn1, seedsfn1, datasmthfn1, cutofffn1, pkcutofffn1] = refine_sig(m, roifn1, bgfn, sigfn1, bgffn, seedsfn1, datasmthfn1, cutofffn1, pkcutofffn1, Puse.p, Puse.options); %% final clean seeds %% sz = Params.neuron_size; [roifn, sigfn, seedsfn, datasmthfn, cutofffn, pkcutofffn] = final_seeds_select(m, roifn1, sigfn1, seedsfn1, datasmthfn1, cutofffn1, pkcutofffn1, sz, imaxy); %% final trace clean %% tflag = 2; sigfn = trace_clean(sigfn, Fsi_new, tflag); %% final refine sig %% [sigfn, spkfn] = pure_refine_sig(sigfn, Puse.options); %% final clean outputs %% sigfn = max(roifn, [], 1)' .* sigfn; roifn = roifn ./ max(roifn, [], 1); % dff = compute_dff(sigfn, bgfn, bgffn, seedsfn); %%% estimate df/f %%% imcur = normalize(imaxy1); imref = normalize(imaxy); [img, sx, sy] = logdemons_unit(imref, imcur); for ii = 1: length(sx) ibmean = iminterpolate(ibmean, sx{ii}, sy{ii}); end x = (imx1 - imn1) * (imx2 - imn2) + imn1; roifnt = roifn; roifnt = roifnt ./ sum(roifnt, 1); bguse1 = ibmean(:)' * roifnt; bguse2 = min(sigfn, [], 2) * x; bguse = bguse1(:) * (imx1 - imn1) + bguse2(:); dff = double(full((sigfn - min(sigfn, [], 2)) * x ./ bguse)); %% save data %% stype = parse_type(class(m.reg(1, 1, 1))); nsize = pixh * pixw * nf * stype; %%% size of single %%% nbatch = batch_compute(nsize); ebatch = ceil(nf / nbatch); idbatch = [1: ebatch: nf, nf + 1]; nbatch = length(idbatch) - 1; imax = zeros(pixh, pixw); for j = 1: nbatch tmp = m.reg(1: pixh, 1: pixw, idbatch(j): idbatch(j + 1) - 1); imax = max(cat(3, max(tmp, [], 3), imax), [], 3); end file_name_to_save = [path_name, file_base{i}, '_data_processed.mat']; if exist(file_name_to_save, 'file') % if ismc % load(file_name_to_save, 'raw_score', 'corr_score') % end delete(file_name_to_save) end if ismc try load(m.Properties.Source, 'raw_score', 'corr_score') save(file_name_to_save, 'roifn', 'sigfn', 'dff', 'seedsfn', 'spkfn', 'bgfn', 'bgffn', 'imax', 'pixh', 'pixw', 'corr_score', 'raw_score', 'Params', '-v7.3'); catch save(file_name_to_save, 'roifn', 'sigfn', 'dff', 'seedsfn', 'spkfn', 'bgfn', 'bgffn', 'imax', 'pixh', 'pixw', 'Params', '-v7.3'); end else save(file_name_to_save, 'roifn', 'sigfn', 'dff', 'seedsfn', 'spkfn', 'bgfn', 'bgffn', 'imax', 'pixh', 'pixw', 'Params', '-v7.3'); end save(file_name_to_save, 'imaxn', 'imaxy', 'imeanf', '-append'); time1 = toc(hpipe); disp(['Done all, total time: ', num2str(time1), ' seconds']) else filename_raw = [path_name, file_base{i}, '_frame_all.mat']; filename_reg = [path_name, file_base{i}, '_reg.mat']; file_name_to_save = filecur; time1 = toc(hpipe); disp(['Done all, total time: ', num2str(time1), ' seconds']) end end end function min1pipe_init % parse path, and install cvx if not % Jinghao Lu, 11/10/2017 %%% prepare main folder %%% pathname = mfilename('fullpath'); mns = mfilename; lname = length(mns); pathtop1 = pathname(1: end - lname); %%% check if on path %%% pathCell = regexp(path, pathsep, 'split'); if ispc % Windows is not case-sensitive onPath = any(strcmpi(pathtop1(1: end - 1), pathCell)); %%% get rid of filesep %%% else onPath = any(strcmp(pathtop1(1: end - 1), pathCell)); end %%% set path and setup cvx if not on path %%% cvx_dir = [pathtop1, 'utilities']; pathcvx = [cvx_dir, filesep, 'cvx', filesep, 'cvx_setup.m']; if ~onPath pathall = genpath(pathtop1); addpath(pathall) if ~exist([cvx_dir, filesep, 'cvx'], 'dir') if ispc cvxl = 'http://web.cvxr.com/cvx/cvx-w64.zip'; elseif isunix cvxl = 'http://web.cvxr.com/cvx/cvx-a64.zip'; elseif ismac cvxl = 'http://web.cvxr.com/cvx/cvx-maci64.zip'; end disp('Downloading CVX'); unzip(cvxl, cvx_dir); end end if ~exist(fullfile(fileparts(prefdir), 'cvx_prefs.mat'), 'file') run(pathcvx); end end