% MATLAB function for calculating the relative contact area after the algorithm in algorithm.txt has been conducted.
% Function originally developed by Arttu Miettinen, Ph.D., modifications done by Tuomas Sormunen, M.Sc.

function plot_contact_area()
    addpath('c:\mytemp\matlab');
    format long
    close all hidden
    samples(1).name = '20171110_160s_b1_1081im';
    samples(1).title = '20171110 160s b1 1081im';
    samples(1).width = 1000;
    samples(1).height = 1025;
    samples(1).pixel_size = 0.065; % µm


    for i = 1:numel(samples)
        
        [results(i).separation, results(i).rac, results(i).thickness,...
            results(i).cumulative_thickness, results(i).limits]...
            = process(samples(i).name, samples(i).title, samples(i).width, samples(i).height, samples(i).pixel_size);
        
    end
    
    figure;
    clf;
    hold on;
    for i = 1:numel(samples)
        plot(results(i).separation, results(i).rac);
    end  
    xlabel('Separation [µm]');
    ylabel('Relative contact area');
    ylim([0, 1]);
    legend({samples(:).title}, 'location', 'southeast');
    topdf(gcf, 'separation_rac.pdf');
    
    
    
    hold off
    figure;
    hold on
    for i = 1:numel(samples)
        plot(results(i).limits, results(i).thickness);
    end
    xlabel('Thickness [µm]');
    ylabel('Relative occurrence');
    legend({samples(:).title}, 'location', 'northwest');
    topdf(gcf, 'thickness.pdf');
    
    
    
    hold off
    figure;
    hold on
    for i = 1:numel(samples)
        plot(results(i).limits, results(i).cumulative_thickness);
    end
    xlabel('Cumulative thickness');
    ylabel('Percentage');
    legend({samples(:).title}, 'location', 'northwest');
    topdf(gcf, 'cumulative_thickness.pdf');
    
    
end

function [separation, relative_contact_area, h, cumulative_h, limits] ...
    = process(name, sample_title, width, height, pixel_size)

    mask = readrawstackd([name '_contact_mask_' num2str(width) 'x' num2str(height) 'x1.raw'], 'uint8');
    opening_thickness = readrawstackd([name '_opening_thickness_' num2str(width) 'x' num2str(height) 'x1.raw'], 'float32');
    thickness = readrawstackd([name '_thickness_' num2str(width) 'x' num2str(height) 'x1.raw'], 'float32');

    mask = mask ~= 0;
    mask = single(mask);
    
    total_area = sum(mask(:) ~= 0);
    
    separation = 1:max(opening_thickness(:));
    contact_area = zeros(size(separation));
    for i = 1:numel(separation)
        contact_area(i) = sum(opening_thickness(:) < separation(i)) - sum(mask(:) == 0);
    end

    relative_contact_area = contact_area ./ total_area;
    separation = separation * pixel_size;
    
    
    
    
    centroid = input('Enter centre of mass of contact mask in the form [x,y] : ');
    
    
    total_area
    %pixel_size
    %abs_total_area = total_area*(pixel_size*10^-6)^2
    rel_co_area = relative_contact_area(1)
    %abs_relative_area = relative_contact_area(1)*total_area*(pixel_size*10^-6)^2
    
    

    
    op_thick = zeros(1, total_area);
    thick = zeros(1, total_area);
    distance = zeros(1, total_area);
    
    index = 1;
    for i = 1:width
        for j = 1:height
            if mask(j,i) == 1
                op_thick(index) = opening_thickness(j,i);
                thick(index) = thickness(j,i);
                distance(index) = sqrt((i-centroid(1))^2+(j-centroid(2))^2);
                index = index + 1;
            end
        end
    end
    

    
    %figure
    %histogram(op_thick)
    %figure

    limits = min(thick(:)):max(thick(:));
    [h,bins]=histcounts(thick,length(limits));
    mean_thickness = mean(thick)
    %plot(limits,h/length(thick))
    h = h/length(thick);
    cumulative_h = zeros(1,length(h));
    cumulative_h(1) = h(1);
    for i=2:length(h)
        cumulative_h(i) = cumulative_h(i-1) + h(i);
    end
    
    r1 = sum((thick-mean(thick)).*(op_thick-mean(op_thick)))/sqrt(sum((thick-mean(thick)).^2).*sum((op_thick-mean(op_thick)).^2))
    r2 = sum((distance-mean(distance)).*(op_thick-mean(op_thick)))/sqrt(sum((distance-mean(distance)).^2).*sum((op_thick-mean(op_thick)).^2))
    r3 = sum((distance-mean(distance)).*(thick-mean(thick)))/sqrt(sum((distance-mean(distance)).^2).*sum((thick-mean(thick)).^2))
    
    
    %scatter(thick, op_thick,'x')
    %r1 = corr2(thick, op_thick)
    %r2 = corr2(distance, op_thick)
    %r3 = corr2(distance, thick)
    
    %figure
    %plot(limits,cumulative_h/length(thick));
    %cumulative_h = cumulative_h/length(thick);
end


