%% configurations
%batch number, for when multiple batches exist within a brain
batch_num=1;
startingsliceidx=1; %Change this each batch, SEQUENTIAL for whole brain

%change this if multiple batches are processed together, so that bardensr
%uses the same basecalling threshold.
use_predefined_thresh=0;
% may need to adjust hyb thresholds.
hybthresh=[50 10 40 10];
hybbgn=[100 60 40 50]; %second channel is higher to remove bleedthrough from the txred channel
% if swapped Alex488 with 594 for hyb channel, then set hyb_reg_ch=1 
hyb_reg_ch=3;

%codebooks
codebook_name='codebookM1all.mat';
codebookhyb_name='codebookhyb.mat';

%bc calling
is_barcoded = 0; %0 if not barcoded, 1 if barcoded
rolthresh = [30 30 30 30];
count_thresh=3;
err_corr_thresh=2;% within a cell body, how to collapse barcodes



% path should not be necessary if BARseq_envs was installed properly
%javaaddpath('C:\barseq_env\FIJI_jar\mij.jar')
%% Logging on
t=datetime('today','Format','yyyyMMdd');
diary([char(t),'_a01.log']);

%% check inputs
if isfile(codebook_name)&&isfile(codebookhyb_name)
    codebook_name=['../',codebook_name];
    codebookhyb_name=['../',codebookhyb_name];
else
    error('Codebooks do not exist.\n')
end

try 
    Miji(false)
    MIJ.exit
    fprintf('MIJ installed correctly.\n')
catch ME
    javaaddpath('C:\barseq_env\FIJI_jar\mij.jar')
    try 
        Miji(false)
        MIJ.exit
    catch ME2
        error('MIJ was not set up and not found in path. Were all env files set up properly?\n')
    end
    warning('MIJ was not set up. Successfully added MIJ to path.\n')
end



%%
cd processed
clc
load('40xto10x.mat','tform40xto10x');
batches=ones(1,numel(tform40xto10x))*batch_num; %Change this based on each batch number! 
save('batches.mat','batches');

bardensr_cmdout=run_bardensr(use_predefined_thresh);
%%
import_bardensr_results(batches, codebook_name);
%% call hyb rolonies
load(codebookhyb_name,'codebookhyb');

basecall_hyb(hybthresh,hybbgn,codebookhyb);

slice=35;
check_hyb_call(slice,codebookhyb)
exportgraphics(gcf,fullfile('..',['FOV',num2str(slice),'_hybcall.jpg']));

%

%% Segmentation with cellpose

% run cellpose in python, dilate each cell by 3 pixels
cellpose_cmdout=run_cellpose();
dilation_radius=3;
import_cellpose(dilation_radius);


%% assign rolonies to cells, transform all coordiates to 10x
assign_gene_rolonies();

rolonies_to_10x();
%
use_mock=1;
calculate_depth(use_mock);



organize_processed_data(startingsliceidx);

%% find overlapping cells and remove them, also set neurons.barcoded if with barcodes.

data_fname=dir('alldata*.mat');
data_fname=sort_nat({data_fname.name});
data_fname=data_fname{end};

boxsize=5; % box size for identifying overlaps
pixelsize=6.5/10; % 10x pixel size
filter_overlapping_neurons(data_fname,boxsize,pixelsize);

%% basecall and assign barcodes, if barcoded
if is_barcoded==1
    gaussrad=0;
    basecall_barcodes(rolthresh, gaussrad); %change how barcodes are "called" - but not for all of them just those in cell
    % assign barcodes rolonies to neurons
    assign_bc2cell(); %%This is the code I probably have to change
    %transform barcodes to 10x coordinates
    bc_to_10x();
    % organize and error correct barcode rolonies
    data_fname=dir('alldata*.mat');
    data_fname=sort_nat({data_fname.name});
    data_fname=data_fname{end};
    mismatch_thresh=1; % allow this mismatch when matching barcodes
    organize_bc_data(count_thresh, err_corr_thresh,data_fname,mismatch_thresh);
    % basecall somas and add soma bc data
    mem=memory;
    thread_num=floor(mem.MemAvailableAllArrays/2^30/8);% need ~8GB ram per thread for basecalling somas.
    basecall_somas(thread_num);
    add_somabc();
    % QC soma barcodes
    filter_somabc('filt_neurons-bc-somas.mat', ...
        'complexity',-0.9, ...
        'score',0.85, ...
        'signal',150 ...
        );
end
%% compress RGB to jpeg, and intermediate image files
convert_RGB_jpg
cleanup_files


%% logging off
diary off