cell_downsample_Hessian_segmentation.m

function cell_downsample_Hessian_segmentation(props)

% =========================================================================
img_fname = props.img_fname;
case_string = props.case_string;
optionFigure = props.optionFigure;
window_range = props.window_range;
downsample_factor = props.downsample_factor;
con = props.connectivity;

peak_min = props.peak_min;
basic_t = props.basic_t;
area_min = props.area_min;
area_max = props.area_max;
% =========================================================================

info = imfinfo(img_fname);
num_images = numel(info);

x1 = window_range(1);
x2 = window_range(2);
y1 = window_range(3);
y2 = window_range(4);

downsize_halfWin = max(floor(downsample_factor/2),1);
ave_filter = fspecial('average', [downsample_factor downsample_factor]);

case_dir = fullfile('cases',case_string);
mkdir(case_dir);
rmdir(case_dir,'s');
mkdir(case_dir);

out_dir_images = fullfile(case_dir,'images');
mkdir(out_dir_images);

out_dir_segmentation = fullfile(case_dir,'segmentation');
mkdir(out_dir_segmentation);

% if optionFigure   
%     
%     out_dir_frames = fullfile(case_dir,'frames');
%     mkdir(out_dir_frames);
%     
%     myVideo = VideoWriter(fullfile(out_dir_frames,sprintf('undersampled_%.2d.avi',downsample_factor)));
%     myVideo.FrameRate = 5;  % Default 30
%     myVideo.Quality = 75;    % Default 75
%     open(myVideo)
%     figure('color','w','position',[0 200 1800 1000])
% end
% 

parfor current_frame = 1:num_images
    
    img_raw = imread(img_fname, current_frame);
    img_cropped = medfilt2(img_raw(y1:y2,x1:x2));
    
    Y = prctile(img_cropped(:),99.9);
    img_adjusted = uint8(double(img_cropped)/double(Y) * 0.9 * double(intmax('uint8')));
    
    img_aveFiltered = uint8(filter2(ave_filter, img_adjusted));
    s = size(img_aveFiltered);
    
    n_x = floor(s(2)/downsample_factor);
    n_y = floor(s(1)/downsample_factor);
    img_aveFiltered = img_aveFiltered(1:n_y*downsample_factor,1:n_x*downsample_factor);
    img_downSampled = img_aveFiltered(downsize_halfWin:downsample_factor:end-downsize_halfWin,downsize_halfWin:downsample_factor:end-downsize_halfWin);
    
    [lambda1,lambda2] = Hessian_2D(img_downSampled);
    He_mask = (lambda2<0)&(lambda1<0);
    He_mask([1:3,end-3:end],:) = 0;
    He_mask(:,[1:3,end-3:end]) = 0;
    cen_He = cell2mat(struct2cell(regionprops(He_mask,'Centroid'))');
    N_mask = length(cen_He);
    
    img_mean = ceil(mean(img_downSampled(:)));
    BW = imbinarize(img_downSampled, img_mean/255);
    cen_BW = cell2mat(struct2cell(regionprops(BW,'Centroid'))');
    [L, N_BW] = bwlabel(BW);
    for r = 1:N_mask
        [dis_min, true_num] = finding_true_num_2D(cen_BW,cen_He(r,:));
        L(L==true_num) = N_BW + 1;
    end
    BW_new = (L == N_BW + 1);
    L_new = bwlabel(BW_new,8);
    
    img_HessianFiltered = img_downSampled;
    img_HessianFiltered(~BW) = 0;
    
    core = cell_core_segmentation(img_HessianFiltered,peak_min,con);
    
    seg_raw = repelem(core,downsample_factor,downsample_factor);
    seg = cell_refinement(seg_raw,img_aveFiltered,basic_t,area_min,area_max);
    parsave(fullfile(out_dir_segmentation,sprintf('segmentation_%.4d',current_frame)),seg)
    
    img_adjusted = img_adjusted(1:n_y*downsample_factor,1:n_x*downsample_factor);
    filename = fullfile(out_dir_images,sprintf('image_%.4d.tif',current_frame));
    imwrite(img_adjusted,filename);   %for saving preprocessed images    
    
%     if optionFigure
%         
%         subplot(361)
%         imshow(img_cropped)
%         title('Raw image')
%         
%         subplot(362)
%         imshow(img_adjusted)
%         title('Intensity adjusted')
%         
%         subplot(363)
%         imshow(img_aveFiltered)
%         title('Moving averaging filtered')
%         
%         subplot(367)
%         imshow(img_downSampled)
%         title(sprintf('Downsampled by factor of %.1d',downsample_factor))
%         
%         subplot(368)
%         imagesc(lambda1)
%         caxis([-10 10])
%         axis equal off
%         title('Hessian \lambda_1')
%         
%         subplot(369)
%         imagesc(lambda2)
%         axis equal off
%         caxis([-10 10])
%         title('Hessian \lambda_2')
%         
%         subplot(3,6,13)
%         imshow(img_HessianFiltered)
%         title('Hessian filtered')
%         
%         subplot(3,6,14)
%         imagesc(core)
%         axis equal off
%         title('Segmented cell core')
%         
%         subplot(3,6,15)
%         imagesc(seg)
%         axis equal off
%         title('Refined cell profile')
%         
%         subplot(3,6,[4 5 6 10 11 12 16 17 18])
%         imshow(img_adjusted); hold on;
%         contour(seg~=0,1,'b','linewidth',2); hold on;
%         %     plot(cell_pos(:,1),cell_pos(:,2),'rX','markersize',5,'linewidth',2,'linesmoothing','on'); hold on;
%         %     ellipse(cell_paras(:,1),cell_paras(:,2),cell_paras(:,3),cell_pos(:,1),cell_pos(:,2),'r',50); hold on;
%         title(sprintf('Segmented cell from frame = %.3d',current_frame))
%         
%         %         print(fullfile(out_dir_frames,sprintf('segmentation_frame=%.4d.tiff',current_frame)),'-dtiff','-r100')
%         set(gca,'nextplot','replacechildren');
%         frame = getframe(gcf);
%         writeVideo(myVideo,  frame)
%     end
end
% if optionFigure
%     close(myVideo)
% end