function [cellid1,seq1,seqC1,sig1,score1]=mmbasecallbccells(filename,nucprofile,c)
%using segmented images, call barcodes from cells based on INTENSITY. If a nucprofile is
%provided, then reduce nuclear background. nucprofile should be a 1x4
%vector that gives the relative strengths of nuc background signals in the
%four channels.
if ~exist('nucprofile','var')
    nucprofile=[];
end

if ~exist('c','var')
    c=0.8;
end

%% Get folders and unzip aligned bc images, cell mask

 img_folder='F:\YW080_30_33\processed';
% analysis_folder=cd(img_folder);
% output_folder=fullfile(analysis_folder,'figs','bc_cell_crops');
% mkdir(output_folder);
% folders=get_folders();

% barseq_filename = 'alignedbcseq01.tif'

    bcseq_filename={};
    %extract all bc images
    for n=1:15
        bcseq_filename{n}=['alignedfixedregn2vbcseq',num2str(n,'%.2u'),'.tif'];
        extract_aligned_img(filename,bcseq_filename{n});
    end

    %extract cell mask for this FOV
    extract_aligned_img(filename,'cellmask.mat');
    cd([img_folder,'\',filename,'\aligned'])
    load('cellmask.mat','maski');



%% Now for find each cell in FOV
cellid1=unique(maski(:));
cellid1=cellid1(cellid1~=0);

% filename = [img_folder,'\',folders{i},'\aligned'];

if ~isempty(cellid1)
    files=dir('*bcseq*.tif');
    files=sort_nat({files.name});
    info=imfinfo(files{1});
    
    im=zeros(info(1).Height,info(1).Width,4);
    for i=1:length(files)
        for n=1:4
            im(:,:,n,i)=imread(files{i},n);
        end
    end
    im=reshape(im,size(im,1)*size(im,2),4,[]);
    
    mask_linear=reshape(maski,[],1);
    
    %read signal, find mean values in each channel for entire cell mask
    %in all cell IDs
    sig1=zeros(size(cellid1,1),size(im,2),size(im,3));
    for i=1:size(sig1,1)
        sig1(i,:,:)=mean(im(mask_linear==cellid1(i),:,:));
    end
    
    %if nucprofile is provided, 
    if ~isempty(nucprofile)
        [minsig,I]=min(sig1,[],2);
        nucsig=repmat(minsig./reshape(nucprofile(I),[],1,size(minsig,3)),1,4,1).*repmat(nucprofile,size(minsig,1),1,size(minsig,3));
        sig1=max(sig1-nucsig*c,0);
    end
    
    %call sequences
    [score1,seq1]=max(sig1,[],2);
    seq1=reshape(seq1,size(seq1,1),[]);
    seqC1=char(seq1);
    seqC1(seq1==1)='G';
    seqC1(seq1==2)='T';
    seqC1(seq1==3)='A';
    seqC1(seq1==4)='C';
    
    score1=score1./sqrt(sum(sig1.^2,2));
    score1=reshape(score1,size(score1,1),[]);
    score1(isnan(score1))=0.5;
else
    sig1=[];
    seq1=[];
    seqC1=char([]);
    cellid1=[];
    score1=[];
end


end

%% functions
% 
function [folders,pos,xxx,yyy]=get_folders()
    %get all folder names and corresponding slice names
    folders=dir('MAX*');
    folders=sort_nat({folders.name});
    pos={};
    xxx=[];
    yyy=[];
    for n=1:numel(folders)
        i1=regexp(folders{n},'_');
        pos{n}=folders{n}(i1(1)+1:i1(2)-1);
        xxx(n)=str2double(folders{n}(i1(2)+1:i1(3)-1));
        yyy(n)=str2double(folders{n}(i1(3)+1:end));
    end
    xxx=xxx+1; %xxx and yyy +1
    yyy=yyy+1; %xxx and yyy +1
    
end

function extract_aligned_img(folder,filename)
warning('off', 'MATLAB:MKDIR:DirectoryExists');

zip_name=fullfile(folder,'aligned.zip');
file_name=fullfile(folder,'aligned',filename);
output_file_name=file_name;

if ~exist(output_file_name,'file')
    mkdir(fullfile(folder,'aligned'));
    unzip_file(zip_name,file_name,output_file_name)
end
warning('on', 'MATLAB:MKDIR:DirectoryExists');

end


function unzip_file(zip_name,file_name,output_file_name)
%windows file paths doesn't work in java, switch to '/'

file_name(file_name=='\')='/';
zip_name(zip_name=='\')='/';
output_file_name(output_file_name=='\')='/';

zipFilename = zip_name;
zipJavaFile  = java.io.File(zipFilename);
zipFile = org.apache.tools.zip.ZipFile(zipJavaFile);
% entries = zipFile.getEntries; % to get all entries then iterate
entry = zipFile.getEntry(file_name); % empty if not found
if(~isempty(entry))
    inputstream = zipFile.getInputStream(entry);
    outJavaFile = java.io.File(output_file_name);
    outStream = java.io.FileOutputStream(outJavaFile);
    copier = com.mathworks.mlwidgets.io.InterruptibleStreamCopier.getInterruptibleStreamCopier;
    copier.copyStream(inputstream, outStream);
    outStream.close();
end
zipFile.close();
end

function crop_im=make_crop_RGB(seqfilename,posx,posy,sz)
max_perct=99.9;

%generate RGB crops of cells for sequencing images
im=double(imread(seqfilename));
for n=2:4
    im(:,:,n)=double(imread(seqfilename,n));
end
%pad image in case it's at the edgedom_bc
im=padarray(im,[sz,sz],'both');%dom_bc
posx=posx+sz;
posy=posy+sz;

cmap=[0 1 1;
    1 1 0;
    1 0 1;
    0 0 1];

im_rgb=reshape(reshape(im,[],4)*double(cmap),size(im,1),size(im,2),[]);

crop_im=im_rgb(round(posy-sz/2):round(posy-sz/2)+sz-1,round(posx-sz/2):round(posx-sz/2)+sz-1,:);
max_thresh=sort(crop_im(:),'ascend');
max_thresh=max_thresh(round(numel(crop_im)*max_perct/100));
crop_im=min(crop_im/max_thresh,1);
end
