cell_parameterization.m

function cell_parameterization(props)

case_string = props.case_string;
optionDebug = props.optionDebug;

case_dir = fullfile('cases',case_string);

if optionDebug
    debug_dir = fullfile(case_dir,'debug');
    mkdir(debug_dir)
end

seg_dir = fullfile(case_dir,'segmentation');
seg_list = dir(fullfile(seg_dir,'*.mat'));

GSM_dir = fullfile(case_dir,'GSM');
mkdir(GSM_dir);

particle_dir = fullfile(case_dir,'particles');
mkdir(particle_dir);

image_dir = fullfile(case_dir,'images');
img_list = dir(fullfile(image_dir,'*.tif'));
num_images = length(img_list);

w = 100;

measures = cell(num_images,1);

for current_frame = 1:num_images
    
    img_raw = imread(fullfile(image_dir,img_list(current_frame).name));    
    
    load(fullfile(seg_dir,seg_list(current_frame).name));
    seg = data;
    
    seg_pad = padarray(seg,[w w],'both');
    img_pad = padarray(img_raw,[w w],'both');
    
    cell_N = max(seg(:));
    measures_currentFrame = zeros(cell_N,27);
    
    parfor j = 1:cell_N
        
        BW_j = (seg_pad == j);
        
        img_j = img_pad;
        img_j(~BW_j) = 0;        
       
        cen = round(cell2mat(struct2cell(regionprops(BW_j,'Centroid'))));
        
        if isempty(cen)
            continue
        end
        
        [y,x] = ind2sub(size(BW_j),find(BW_j));
        ww = max(max(abs(y - cen(2))), max(abs(x - cen(1)))) + 25;
      
        x1 = cen(1)-ww;
        x2 = cen(1)+ww;
        y1 = cen(2)-ww;
        y2 = cen(2)+ww;        
        
        img_crop = double(img_j(y1:y2,x1:x2));        
        BW_crop_edge = edge(BW_j(y1:y2,x1:x2),'log');
        BW_crop = imfill(BW_crop_edge,'holes');
        area = cell2mat(struct2cell(regionprops(BW_crop,'area')));
        if length(area) > 1            
            [~,max_ind] = max(area);
            BW_crop = (bwlabel(BW_crop) == max_ind);
        end
        
        %                     figure('position',[0 0 900 900])
        %                     subplot(231)
        %                     imshow(img_raw_crop)
        %                     subplot(234)
        %                     imagesc(img_crop); axis equal off
        %                     subplot(232)
        %                     imshow(BW_crop)
        %                     subplot(233)
        %                     imshow(BW_crop_edge)
        %                     subplot(232)
        %                     imshow(BW_crop)
        %                     subplot(233)
        %                     imshow(BW_crop_edge)
        
        measure = calculate_GSM(img_crop,BW_crop,optionDebug);
        measure(1) = measure(1) + x1 - w;
        measure(2) = measure(2) + y1 - w;
        measure(3) = measure(3) + x1 - w;
        measure(4) = measure(4) + y1 - w;
        if find(inf(isnan(measure) | isinf(measure)))
            continue
        end
        measures_currentFrame(j,:) = measure;
     
%         if optionDebug
%             
%             [y,x] = ind2sub(size(BW_crop_edge),find(BW_crop_edge));
%             
%             subplot(241)
%             imshow(img_crop); hold on;
%             plot(x,y,'.r');
%             title('MFI image')            
%            
%             
%             print(fullfile(debug_dir,sprintf('_%.02d.tif',j)),'-dtiff','-r200')
%         end
        
    end
    ind = measures_currentFrame(:,1) == 0;
    measures_currentFrame(ind,:) = [];
    save(fullfile(GSM_dir,sprintf('%0.4d',current_frame)),'measures_currentFrame')    
    measures{current_frame} = measures_currentFrame;
end

save(fullfile(case_dir,'GSM_all'),'measures')    
measures_all = cell2mat(measures);

cats = {'IWC','Center of mass',...
    'I_{mean}','I_{max}','Area (pixel^2)', 'Perimeter (pixel)',...
    'Major axis (pixel)', 'Minor axis (pixel)','Equiv R (pixel)', 'Orientation (°)', ...
    'Elongation (FFT)','Ecc (chord)','Ecc (ellipse)', 'Elongatedeness', 'Ratio (Rb/Ra)',...
    'Circ (AP)', 'Circ (r_{std})', 'Circ (r_{avg})', 'Circ (r_{max})', ...
    'Ellipticity','Squareness', 'Triangularity',  'Rectangularity', 'Convexity (A)', 'Convexity (P)'};
sub_pos = [0 0 0 0 1 2 3 4 7 8 9 10 13 14 15 16 17 19 20 21 22 25 26 27 28 29 30];
close all;
% figure('color','white','position',[0 0 1800 1000])
% for i = 5:27
%     subplot(5,6,sub_pos(i))
%     measure = measures_all(:,i);  
%     ind = find(isnan(measure) | isinf(measure));
%     measure(ind) = [];    
%     measure_avg = abs(mean(measure));
%     measure_std = var(measure).^0.5;
%     edges = linspace(measure_avg-3*measure_std,measure_avg+3*measure_std,20);
%     histogram(measure,edges,'Normalization','probability')
%     xlim([min(measure) max(measure)])
%     title({cats{i-2},sprintf('RMS/MEAN = %.2f',measure_std/measure_avg)})
%     xlim([measure_avg-3*measure_std measure_avg+3*measure_std])
% end
% print(fullfile(case_dir,'GSM_pdf.tif'),'-dtiff','-r300')


measures_ind = [10,9,12,6,7,16,19,23];
for current_frame = 1:num_images
    
    measures_currentFrame = measures{current_frame};
    N = length(measures_currentFrame);  
    pos = cat(2,measures_currentFrame(:,1:2),zeros(size(measures_currentFrame,1),1));
    paras = measures_currentFrame(:,measures_ind);    
    paras(:,3) = pi/2 - paras(:,3)/180 * pi;
    particles = [];
    particles.pos = pos;
    particles.paras = paras;
    save(fullfile(particle_dir,sprintf('%0.4d',current_frame)),'particles')
    
end

%
%
% particles = [];
% particles.pos = cell_pos;
% particles.paras = cell_paras;
% file_name = fullfile(out_dir_particles,sprintf('particles_%.4d.mat',current_frame));
% parsave(file_name,particles)