%% 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=[50 60 40 30]; %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;
rolthresh = [30 30 30 30];
count_thresh=3;
err_corr_thresh=2;% within a cell body, how to collapse barcodes



load chprofile20x-cypress-20231122.mat
load chshift20x-scopeCypress-20231218.mat
chshift20x=chshift_20x';
 
chshift20x=chshift20x(1:4,1);%in case this contains 5 channels.

%% Logging on
t=datetime('today','Format','yyyyMMdd');
diary([char(t),'_a00.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

%%



%% Preprocess images
clc
numcores=feature('numcores');

p=gcp('nocreate');
if isempty(p)
    parpool("local",round(numcores*1.5));
elseif p.NumWorkers<round(numcores*1.5)
    parpool("local",round(numcores*1.5));
end

%% Process geneseq

organize_geneseq; %move files
cmdout=n2v_processing('n2vprocessing.py'); %n2v denoise

% Load channel profile for 20x

%%
cd processed
%register geneseq files
fname='n2vgene';
ball_radius=6;
rgb_intensity=0.6;
tic
register_seq_images(fname,chprofile20x,ball_radius,chshift20x,rgb_intensity)
t=toc;

% optional, to remake rgb files
%remake_rgb_img(fname,rgb_intensity)
%% 

%start here
% Make codebook for bardensr
make_codebook(codebook_name)
%% 

%% Process hyb images
organize_hyb_files_multi();%move files
%%
cd ..
n2v_processing('n2vprocessing_hyb.py');%n2v hyb images
cd processed
%%
chprofilehyb=[1 0 0 0;0.52 1 0 0;0 0 1 0;0 0 0 1]; % measured on this dataset
chshifthyb=[0 -3;0 0;0 0;-1 4]; % measured on this dataset
save('chprofilehyb.mat','chprofilehyb','chshifthyb');

% register images
%hyb_reg_ch=3; %channel used to align
nuclearch=5; %nuclear channel
radius=100; %radius for bgn subtraction of all channels.
reg_cycle=1; % register to 1st cycle
chradius=30; % bgn radius for seq rolony cycle
rgb_intensity=0.6;

register_hyb_images_multi(hyb_reg_ch,radius, nuclearch,reg_cycle,chradius,chprofilehyb,chshifthyb, rgb_intensity);

%% Stitch images from first cycle
ch_count=5;
reg_ch=4;
overlap=0.23;
x=1;
trim=0.02;
rescale_factor=0.5; % rescaling images by this
stitch_10x_images(ch_count,reg_ch,overlap,x,trim,rescale_factor);

%% generate images to visually check registration quality        
check_stitching('40xto10x.mat')
% stitch the transformed 40x image and overlay with 10x images for visual inspection. This produces 5x smaller images than the original 10x to ease loading, but should still be enough since the registration only need to be roughly correct.
intensity_scaling=3;
checkregistration('geneseq','40xto10x.mat',intensity_scaling);
checkregistration('hyb','40xto10x.mat',intensity_scaling);

% Here check the images in "chekregistration" to see if all the images match within each pos.        

%%
cd ..

% ORganize barcodes. 
if is_barcoded

    organize_bcseq()
    cmdout=n2v_processing('n2vprocessing_bc.py');
    cd processed
    
    %register BC to geneseq
    BC_refch=5;
    gene_refch=5;
    BC_name='n2vbcseq';
    gene_name='n2vgeneseq';
    alignBC2gene(BC_refch,gene_refch,BC_name,gene_name);
    
    % Register BC 40x
    % Load channel profile for 20x
    load chprofile20x-50-30-20-40-20220218.mat
    load chshift20x-20220218
    chshift20x=chshift20x(1:4,:);%in case this contains 5 channels.
    chshift20x=chshift20x+[1,-0.7;1,0;0,0;0,0]; %manually fixed based on current data.
    
    fname='regn2vbc';
    ball_radius=100;
    rgb_intensity=0.6;
    register_seq_images(fname,chprofile20x,ball_radius,chshift20x,rgb_intensity)
    %
    
    %generate checkregistration files
    intensity_scaling = 3;
    checkregistration('bcseq','40xto10x.mat',intensity_scaling);
    
    % Stitch BC image, this is optional, to produce a first-cycle BC
    %stitched image
    ch_count=5;
    reg_ch=4;
    overlap=0.23;
    x=1;
    trim=0.02;
    rescale_factor=0.5; % rescaling images by this
    niestitchgenessingleplane(ch_count,reg_ch,overlap,x,trim,'regn2vbc');
    
    cd ..
end


%% logging off
diary off