calculate_GSM.m

function measures = calculate_GSM(raw_img,BW_img,optionDebug)
%==========================================================================
% Global Shape Measurements:
% Quantifies region properties of binary images
% of arbitrarily-shaped particulates in auto injectors or living cells
%
% References:
% 1. Carosone, F., Cenedese, A., & Querzoli, G. (1995). Recognition of
% partially overlapped particle images using the Kohonen neural network.
% Experiments in Fluids, 19(4), 225-232.
% 2. Rosin, P. L. (2005). Computing global shape measures.
% In Handbook of Pattern Recognition and Computer Vision (pp. 177-196).
% 3. Horej�, J. (2010). Shape Analysis Using Global Shape Measures.
%
% Version 1.0
% Jan.2019. by Tianqi Guo
%==========================================================================

% basics
CoM = cell2mat(struct2cell(regionprops(BW_img,'centroid'))');
area = cell2mat(struct2cell(regionprops(BW_img,'area')));
r_eq = cell2mat(struct2cell(regionprops(BW_img,'EquivDiameter')))/2;
per = cell2mat(struct2cell(regionprops(BW_img,'perimeter')));
circularity_AP = 4*pi*area/per^2;

% intensity based
IWC = cell2mat(struct2cell(regionprops(BW_img,raw_img,'WeightedCentroid')));
I_avg = cell2mat(struct2cell(regionprops(BW_img,raw_img,'MeanIntensity')));
I_max = cell2mat(struct2cell(regionprops(BW_img,raw_img,'MaxIntensity')));
%==========================================================================
% boundary based

boundary = cell2mat(bwboundaries(BW_img));
boundary_x = boundary(:,2);
boundary_y = boundary(:,1);
n_boundary = length(boundary);

% longest chords
[BW_img_rotated,chords_xy] = longest_chord(BW_img,boundary_x,boundary_y);
chord_shrt = abs(chords_xy(1,3)-chords_xy(2,3));
chord_long = abs(chords_xy(1,2)-chords_xy(2,2));
eccentricity_chord = chord_shrt / chord_long;

% complex Fourier transform
c = abs(fft(boundary_x + 1i * boundary_y)/n_boundary);
elongation = c(end)/c(2);
squareness = c(n_boundary-2)/c(2);

% m = 5;
% ind = (n_boundary - 2 - m*4) : 4 : (n_boundary - 2);
% s_3 = sum(c(ind)/c(2))/sum(c(3:end-1)/c(2));

% =========================================================================
% moments based

% raw moments
M00 = raw_moments(BW_img,0,0);
M10 = raw_moments(BW_img,1,0);
M01 = raw_moments(BW_img,0,1);
M11 = raw_moments(BW_img,1,1);
M20 = raw_moments(BW_img,2,0);
M02 = raw_moments(BW_img,0,2);

% central moments
x_0 = M10/M00;
y_0 = M01/M00;
mu00 = M00;
mu11 = M11 - x_0 * M01;
mu20 = M20 - x_0 * M10;
mu02 = M02 - y_0 * M01;

% Elongatedness
elongatedness = sqrt((mu20 - mu02)^2 + 4 * mu11^2)/(mu20 + mu02);

% moment invariant
I1 = (mu20 * mu02 - mu11^2)/mu00^4;

% triangularity and ellipticity
if I1 <= 1/108
    triangularity = 108 * I1;
else
    triangularity = 1/(108*I1);
end
if I1 <= 1/(16*pi^2)
    ellipticity = 16*pi^2 * I1;
else
    ellipticity = 1/(16*pi^2*I1);
end

% based on equivalent ellipse with the same second-moments
ra = cell2mat(struct2cell(regionprops(BW_img,'MajorAxisLength')))/2;
rb = cell2mat(struct2cell(regionprops(BW_img,'MinorAxisLength')))/2;
eccentricity_ellipse = cell2mat(struct2cell(regionprops(BW_img,'Eccentricity')));
orient_degree = cell2mat(struct2cell(regionprops(BW_img,'Orientation')));
orient_radian = orient_degree/180*pi;
aspectRatio = rb/ra;

%==========================================================================
% bounding region based

% smallest bounding convex polygon
Con_img = cell2mat(struct2cell(regionprops(BW_img,'ConvexImage')));
Con_hul = cell2mat(struct2cell(regionprops(BW_img,'ConvexHull')));
Con_per = cell2mat(struct2cell(regionprops(Con_img,'perimeter')));
convexity_a = cell2mat(struct2cell(regionprops(BW_img,'solidity')));
convexity_p = Con_per/per;

% minimum bounding box (rectangle)
img_rotated = imrotate(double(raw_img), - orient_degree, 'crop','bilinear');
BW_rotated_edge = edge(img_rotated,'log');
BW_rotated = imfill(BW_rotated_edge,'holes');
area_rotated = cell2mat(struct2cell(regionprops(BW_rotated,'area')));

s = regionprops(BW_rotated,'BoundingBox');
s = cat(1, s.BoundingBox); %[left bottom width height]
if size(s,1) > 1
    t = [s(:,1) s(:,2) s(:,1)+s(:,3) s(:,2)+s(:,4)];  %[left bottom right top]
    leftBottomCorner = min(t(:,1:2));
    rightTopCorner = max(t(:,3:4));
    MBB_xy = [leftBottomCorner rightTopCorner-leftBottomCorner]; %[left bottom width height]
else
    MBB_xy = s;
end

rectangularity  = sum(area_rotated) / (MBB_xy(3)*MBB_xy(4));

%==========================================================================
% Polar transform based

% Polar transform
nrad = 300;
ntheta = 360;
linlog = 'linear';
shape = 'valid';
[theta, radius, pim] = polartrans(BW_img, nrad, ntheta, CoM(1), CoM(2), linlog, shape);

% r(thetha)

r_theta = zeros(1,ntheta);
for t = 1:ntheta
    %ind = find(pim(:,t)==0,1,'first');
    ind = find(pim(:,t)==1,1,'last');
    if isempty(ind)
        r_theta(t) = 0;
    else
        r_theta(t) = radius(ind(1),t);
    end
end
r_theta = smooth(r_theta(:));
r_avg = mean(r_theta);
r_std = var(r_theta).^0.5;
r_max = max(r_theta);

% circularities
circularity_Ravg = 1 - trapz(theta(1,:),(r_theta - r_avg).^2)/(r_avg^2*trapz(theta(1,:),ones(size(theta(1,:)))));
circularity_Rmax = area/(pi*r_max^2);
circularity_Rstd = 1 - r_std/r_avg;

% =========================================================================
% if orient_degree < 0
%     orient_degree = orient_degree + 90;
% end
measures = [IWC(1), IWC(2), CoM(1), CoM(2), I_avg, I_max...
    area, per, ra, rb, r_eq, orient_degree,...
    elongation, eccentricity_chord, eccentricity_ellipse, elongatedness, aspectRatio, ...
    circularity_AP, circularity_Rstd, circularity_Ravg, circularity_Rmax,  ...
    ellipticity, squareness, triangularity, rectangularity, convexity_a, convexity_p];
cats = {'Basics','IWC','Center of mass','I_mean','I_max','Area (pixel^2)', 'Perimeter (pixel)',...
    'Major axis (pixel)', 'Minor axis (pixel)','Equiv R (pixel)', 'Orientation (°)', ...
    'Eccentricities','Elongation (FFT)','Ecc (chord)','Ecc (ellipse)', 'Elongatedeness', 'Ratio (Rb/Ra)',...
    'Circularities','Circ (AP)', 'Circ (r_std)', 'Circ (r_avg)', 'Circ (r_max)', ...
    'Shape factors','Ellipticity','Squareness', 'Triangularity',  'Rectangularity', 'Convexity (A)', 'Convexity (P)'};
table_data = {};
table_data{1} = [];
table_data{2} = sprintf('(%.1f,%.1f)',IWC);
table_data{3} = sprintf('(%.1f,%.1f)',CoM);
for i = 4 : 11
    table_data{i} = sprintf('%.1f',measures(i+1));
end
table_data{12} = [];
for i = 13 : 17
    table_data{i} = sprintf('%.4f',measures(i));
end
table_data{18} = [];
for i = 19 : 22
    table_data{i} = sprintf('%.4f',measures(i-1));
end
table_data{23} = [];
for i = 24 : 29
    table_data{i} = sprintf('%.4f',measures(i-2));
end
if optionDebug
    file_list = dir(fullfile('debug','*.tif'));
    N = length(file_list) + 1;
    
    close all;
    fig = figure('position',[00 200 1500 900]);
    
    %     subplot(241)
    %     imshow(BW_rotated);  hold on;
    %     rectangle('Position',MBB_xy,'edgecolor','m','linestyle',':','linewidth',2);
    %     title('Minimum bounding rectangle')
    %
    %     subplot(242)
    %     imshow(BW_img_rotated); hold on;
    %     plot(chords_xy(:,1),chords_xy(:,2),'-m','linewidth',2);
    %     plot(chords_xy(:,3),chords_xy(:,4),'-m','linewidth',2)
    %     title('Longest chord')
    
    subplot(243)
    s = size(raw_img);
    imagesc(uint8(raw_img)); axis equal off; hold on;
    h1 = plot(CoM(1),CoM(2),'bO','markersize',5,'linewidth',2);
    h2 = ellipse(rb,ra,pi/2-orient_radian,CoM(1),CoM(2),'r',50); hold on;
    h3 = plot(Con_hul(:,1),Con_hul(:,2),'g-.','linewidth',2);
    h4 = plot([CoM(1) s(2)],[CoM(2) CoM(2)],'m--','linewidth',2); hold on
    plot([CoM(1) s(2)],[CoM(2) s(1)- (CoM(2)+(s(2)-CoM(1))*tan(orient_radian))],'m--','linewidth',2);
    title('Segmented image')
    legend([h1,h2,h3,h4],...
        {'Center of mass','Equivalent eclipse','Smallest enclosing convex','Orientation angle'},...
        'location','southeast')
    xlim([1 s(2)])
    ylim([1 s(1)])
    
    subplot(245)
    contourf(theta,radius,pim,'edgecolor','none'); hold on;
    plot(theta(1,:),r_theta,'r-','linewidth',2)
    xlabel('\theta (Radian)')
    ylabel('r (pixel)')
    title('Polar transform')
    
    subplot(246)
    plot(boundary_x,'.'); hold on;
    plot(boundary_y,'.'); hold on;
    grid on
    title('Boundary coordinates')
    legend('x','y','location','south')
    xlabel('Number of points')
    ylabel('Coordinates (pixel)')
    xlim([1,n_boundary])
    
    subplot(247)
    n_2 = floor(length(c)/2);
    if mod(n_boundary,2) == 0
        n = -n_2:(n_2-1);
    else
        n = -n_2:n_2;
    end
    r_0 = (x_0^2 + y_0^2)^0.5;
    stem(n,fftshift(c),'filled','linewidth',2,'markersize',5)
    xlim([-10 10])
    ylim([0 r_0*1.1])
    grid on;
    xlabel('Descriptor number (n)')
    ylabel('|a_n + i * b_n|')
    title('Complex Fourier transform')
    
    subplot(2,4,[4,8])
    %cnames = {'Measures','Values'};
    t = uitable(fig,'Data',cat(2,cats(:),table_data(:)),...
        'ColumnName',[],'RowName',[],'ColumnWidth',{120,100},...
        'ColumnFormat',{'char','char'});
    t.FontSize = 10;
    plot(3)
    pos = get(subplot(2,4,[4,8]),'position');
    delete(subplot(2,4,[4,8]))
    set(t,'units','normalized')
    set(t,'position',pos)
    %     jScroll = findjobj(t);
    %     jTable  = jScroll.getViewport.getView;
    %     jTable.setAutoResizeMode(jTable.AUTO_RESIZE_SUBSEQUENT_COLUMNS);
    %     drawnow;
    
    %     print(fullfile('debug',sprintf('%0.4d.tif',N)),'-dtiff','-r100')
    
end
end

function M = raw_moments(BW_img,p,q)
s = size(BW_img);
[X,Y] = meshgrid(1:s(2),1:s(1));
integrand = X.^p .* Y.^q .* BW_img;
M = sum(integrand(:));
end

function [BW_img_rotated,chords_xy] = longest_chord(BW_img,boundary_x,boundary_y)
%==========================================================================
% try
s = size(BW_img);
x0 = round(s(2)/2);
y0 = round(s(1)/2);
n_boundary = length(boundary_x);

chord_length_max = -1;
for i = 1:n_boundary-1
    for j = i+1:n_boundary
        flag = 0;
        if boundary_x(i) == boundary_x(j)
            step = sign(round(boundary_y(j))-round(boundary_y(i)));
            for yy = round(boundary_y(i)) : step : round(boundary_y(j))
                if ~BW_img(yy,boundary_x(j))
                    flag = 1;
                    break
                end
            end
        else
            k = (boundary_y(i) - boundary_y(j))/(boundary_x(i) - boundary_x(j));
            step = sign(round(boundary_x(j))-round(boundary_x(i)));
            for xx = round(boundary_x(i)) : step : round(boundary_x(j))
                yy = round(boundary_y(i) + k * (xx - boundary_x(i)));
                if ~BW_img(yy,xx)
                    flag = 1;
                    break
                end
            end
        end
        if flag
            continue
        end
        
        chord_length = ((boundary_x(i) - boundary_x(j))^2 + (boundary_y(i) - boundary_y(j))^2).^0.5;
        if chord_length > chord_length_max
            chord_length_max = chord_length;
            i_max = i;
            j_max = j;
        end
    end
end
chord_max_x = [boundary_x(i_max) boundary_x(j_max)];
chord_max_y = [boundary_y(i_max) boundary_y(j_max)];

if chord_max_x(1) == chord_max_x(2)
    ang = 0;
elseif chord_max_y(1) == chord_max_y(2)
    ang = pi/2;
else
    ang = pi/2 - atan((chord_max_y(2) - chord_max_y(1))/(chord_max_x(2) - chord_max_x(1)));
end
BW_img_rotated = imrotate(BW_img,-ang/pi * 180,'crop');
chord_max_x_R = round((chord_max_x - x0) * cos(ang) - (chord_max_y - y0) * sin(ang) + x0);
chord_max_y_R = round((chord_max_x - x0) * sin(ang) + (chord_max_y - y0) * cos(ang) + y0);
chord_max_x_R_mid = round(mean(chord_max_x_R));
chord_max_ppd_max = -1;

for yy = min(chord_max_y_R):max(chord_max_y_R)
    
    left = BW_img_rotated(yy,1:chord_max_x_R_mid);
    x1 = find(left==0,1,'last');
    right = BW_img_rotated(yy,chord_max_x_R_mid:s(2));
    x2 = find(right==0,1,'first') + chord_max_x_R_mid;
    if x2 - x1 > chord_max_ppd_max
        chord_max_ppd_max = x2 - x1;
        chord_max_ppd_max_x1 = x1;
        chord_max_ppd_max_x2 = x2;
        chord_max_ppd_max_yy = yy;
    end
end

chords_xy = cat(2,chord_max_x_R(:),chord_max_y_R(:),...
    [chord_max_ppd_max_x2;chord_max_ppd_max_x1],[chord_max_ppd_max_yy;chord_max_ppd_max_yy]);
% catch
%     keyboard
% end
end