%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Treatment of light reflections in 3D PIV systems                %%%
%%% A. Grille Guerra, L. Porcar Galan, A. Sciacchitano & F. Scarano %%%
%%% Delft University of Technology                                  %%%
%%% Measurement Science and Technology, 2025                        %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Image preprocessing and masking routine

% Applied to images from:
% Jux, C., Sciacchitano, A., & Scarano, F. (2020). 
% Flow pressure evaluation on generic surfaces by robotic volumetric PTV. 
% Measurement Science and Technology, 31(10). 
% https://doi.org/10.1088/1361-6501/ab8f46

clear;
close all;
clc;

%% Preliminary 
% Plot parameters
width = 6;     % Width in inches 8
height = 5;    % Height in inches 4
alw = 1;       % AxesLineWidth 
fsz = 20;      % Fontsize 
lw = 2;        % LineWidth 
fsl = 24;      % Font size Label 
msz = 8;       % MarkerSize

%% load data
num_images = 20;
unzip('sphere images.zip');
data_path = '.\sphere images';

% Cut-off frequencies
fc = 0.01; % temporal cut-off, relative to the Nyquist frequency
kc = 0.01; % spatial cut-off, in 1/px
N_filt = 3; % filter order for the Butterworth temporal filter

frame = imread([data_path '\B0001.tif']);
[J,I] = size(frame);
Sequence = zeros(J,I,num_images);

for i=1:num_images
    frame = imread([data_path '\B' num2str(i,'%04d') '.tif']);
    Sequence(:,:,i) = frame;
end
Sequence = uint16(Sequence);

%split images retaining first camera only
Sequence = Sequence(1:end/4,:,:);

%% Temporal minimum subtraction
Min_subtracted = uint16(Sequence-min(Sequence,[],3));

%% High-pass filtering (spatial)
sigma = 2.5; %size=2*ceil(2*sigma)+1
for i=1:num_images
    High_pass(:,:,i) = Sequence(:,:,i)-imgaussfilt(Sequence(:,:,i),sigma);
end

%% Butterworth + high-pass spatial filter (in Fourier space)

%butterworth
[J,I,N] = size(Sequence);
Butterworth = zeros(J,I,N);
[b,a] = butter(N_filt,fc,'high');
if mod(N,2) %odd case
    [h,~] = freqz(b,a,(N+1)/2);
h_both = [flipud(conj(h)); h(2:end)];
else %even case
    [h,~] = freqz(b,a,N/2+1);
h_both = [flipud(conj(h)); h(2:end-1)];
end

for j=1:J
    for i=1:I
        Raw_time_history = squeeze(Sequence(j,i,:));
        Raw_time_history = double(Raw_time_history);
        Y = fftshift(fft(Raw_time_history));
        Filt_time_history = ifft(ifftshift(Y.*abs(h_both)),'symmetric');
        Butterworth(j,i,:) = Filt_time_history;
    end
end
Butterworth = uint16(Butterworth);

%high-pass spatial filtering
for i=1:num_images
    ft_image = fftshift(fft2(Butterworth(:,:,i)));
    x = linspace(-0.5,0.5-1/I,I);
    y = linspace(-0.5,0.5-1/J,J);
    [X,Y] = meshgrid(x,y);
    filter = 1 - exp(-X.^2/(2*kc^2)).* exp(-Y.^2/(2*kc^2));
    G = ft_image.*filter;
    ImPrePro(:,:,i) = uint16(ifft2(ifftshift(G),'symmetric'));
end

%% SSIM masking

% Masking parameters
thresholdSSIM = 0.5;
thresholdCntArea = 50;
for i=1:num_images
    % Obtain SSIM map
    [ssimval,ssimmap] = ssim(Sequence(:,:,i),ImPrePro(:,:,i));
    SSIM(:,:,i) = ssimmap;
    % Get masked image from SSIM coefficient plot by thresholding
    ssimmap(ssimmap>thresholdSSIM) = 1;
    ssimmap(ssimmap<=thresholdSSIM) = 0;
    % Retrieve contours from Structural similarity (SSIM) plot
    % between output and input images
    contours = ssimmap == 0;
    [L, num] = bwlabel(contours, 8); 
    % Count of pixels in each blob (area of each blob)
    counts = sum(bsxfun(@eq,L(:),1:num)); 
    % Find contours biggest than a set threshold area
    [~,ind] = find(counts > thresholdCntArea);
    % Build mask with largest contours
    num_indxs = size(ind,2);
    biggestContours = ones(J, I);
    for j = 1:num_indxs
        K = L;
        out = (K==ind(j));
        if j == 1
            biggestContours = out & biggestContours;
        else
            biggestContours = out + biggestContours;
        end
    end
    % Create mask with biggest contours found
    mask = uint16(~biggestContours);
    Final(:,:,i) = ImPrePro(:,:,i).*mask;
end


%% Plot
frame = 10;
axis_lim = 500;

%circle
r = 120;
x = 450;
y = 240;
th = 0:pi/50:2*pi;
xunit = r * cos(th) + x;
yunit = r * sin(th) + y;
xunit(xunit>568) = NaN;

%sting
x_s = [567 640];
y_top = [220 220];
y_low = [260 260];

%mask
mask = edge(biggestContours);
[row,col] = ind2sub([J,I],find(mask));

% Figure 1: raw, time minimum and high-pass filtered

figure

ax1 = subplot(1,3,1);
imshow(Sequence(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
clim([0 axis_lim])
colormap(flipud(gray(10)))

ax2 = subplot(1,3,2);
imshow(Min_subtracted(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
clim([0 axis_lim])
colormap(flipud(gray(10)))

ax3 = subplot(1,3,3);
imshow(High_pass(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
clim([0 axis_lim])
colormap(flipud(gray(10)))

pos = get(gcf, 'Position');
set(gcf, 'Position', [pos(1)-400 pos(2)-300 width*260, height*90]); %<- Set size
ax1.Position = [0 0 0.33 1];
ax2.Position = [0.335 0 0.33 1];
ax3.Position = [0.67 0 0.33 1];

% Figure 2: pre-processed, SSIM and masked image

figure

ax1 = subplot(1,3,1);
imshow(ImPrePro(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
clim([0 axis_lim])
colormap(flipud(gray(10)))

ax2 = subplot(1,3,2);
imshow(SSIM(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
plot(col,row,'r.','LineWidth',2)
clim([0 1])
colormap(ax2,gray(10))

ax3 = subplot(1,3,3);
imshow(Final(:,:,frame))
hold on
plot(xunit, yunit,'--','LineWidth',3,'Color','g');
plot(x_s, y_top,'--','LineWidth',3,'Color','g');
plot(x_s, y_low,'--','LineWidth',3,'Color','g');
plot(col,row,'r.','LineWidth',2)
clim([0 axis_lim])
colormap(ax3,flipud(gray(10)))

pos = get(gcf, 'Position');
set(gcf, 'Position', [pos(1)-400 pos(2)-300 width*260, height*90]); %<- Set size
ax1.Position = [0 0 0.33 1];
ax2.Position = [0.335 0 0.33 1];
ax3.Position = [0.67 0 0.33 1];