How to implement a multitemporal cloud detection algorithm using Sentinel-2 time series

Hello,

I am having some troubles trying to implement a multitemporal cloud detection algorithm (MTCD) for Sentinel 2 time series images (I have a time series composed of 36 images). I am following the algorithm poposed in a paper that explains the following steps:

1) Make a background image;

2)Create a difference image between the target (cloudy) image and the background;

3)Apply a k-means clustering algorithm to divide the image into 10 clusters;

4)Analyze the clusters by applying to each cluster 3 thresholds:

-a = the norm (intensity) of the difference reflectance image over the visible bands (R,G,B);

-b = the mean of the difference reflectance image over the visible bands (R,G,B);

-y = the norm of the reflectance image over the visible bands (R,G,B).

5)A cluster is then classified as cloudy if the following are satisfied: a>=0.04, b>=0, y>=0.175

My problems start in the fourth step; I am not able to compute the thresholds. I am adding the code I wrote below.

Please can somebody help me with the calculation of the tresholds at a cluster level?

Thank you so much :)

clc;

clear;

images = dir('S2*2017*.tif'); %.tif format of the time series

%-------------------------------------------

% Loading the images bands

%-------------------------------------------

img = single(imread(images(1).name));

X_im_red(:,:,1) = img(:,:,1) ;

%figure,imshow(X_im_red(:,:,1));

X_im_green(:,:,1) = img(:,:,2);

%figure,imshow(X_im_green(:,:,1));

X_im_blue(:,:,1) = img(:,:,3);

%figure,imshow(X_im_blue(:,:,1));

%I still have to load the missing 7 bands

for i=2 :size(images,1)

img = single(imread(images(i).name));

X_im_red(:,:,i) = img(:,:,1);

X_im_green(:,:,i) = img(:,:,2);

X_im_blue(:,:,i) = img(:,:,3);

end

[nl,nc,nb] = size(X_im_blue);

% nl(rown)=1340 nc(colums)=1432 nb(images)=36

%-------------------------------------------

% Background

%-------------------------------------------

X_red = reshape(X_im_red,[nl*nc,nb]);

X_Med(:,1) = quantile(X_red,0.25,2);

X_green = reshape(X_im_green,[nl*nc,nb]);

X_Med(:,2) = quantile(X_green,0.25,2);

X_blue = reshape(X_im_blue,[nl*nc,nb]);

X_Med(:,3) = quantile(X_blue,0.25,2);

X_im_Med = reshape(X_Med,[nl,nc,3]);

figure,imshow(uint8(X_im_Med));

X_im_Med2 = uint8(X_im_Med);

imwrite(X_im_Med2,'background.png','png')

%-------------------------------------------

% Difference Image

%-------------------------------------------

X_im_diff_red = zeros(size(X_im_red));

X_im_diff_green = zeros(size(X_im_green));

X_im_diff_blue = zeros(size(X_im_blue));

for i = 1 : size(X_im_blue,3)

X_im_diff_red(:,:,i) = (X_im_red(:,:,i) - X_im_Med(:,:,1));

X_im_diff_green(:,:,i) = (X_im_green(:,:,i) - X_im_Med(:,:,2));

X_im_diff_blue(:,:,i) = (X_im_blue(:,:,i) - X_im_Med(:,:,3));

end

%I am selecting only the bands belonging to the eighth image (for simplicity)

Xim_red = X_im_diff_red(:,:,8);

Xim_green = X_im_diff_green(:,:,8);

Xim_blue = X_im_diff_blue(:,:,8);

X_dif = cat(3,Xim_red,Xim_green,Xim_blue);

[nl,nc,nb] = size(X_dif); %1340,1432,3

X_dif1 = reshape(X_dif,[nl*nc,nb]);

%-------------------------------------------

% Clustering

%-------------------------------------------

label = kmeans(X_dif1(:),10);

Xlabel = reshape(label,[nl,nc,3]);

color_map = zeros(10,3);

for i=1:10

color_map(i,:) = rand(1,3);

end

rgb = cat(3,X_im_red(:,:,8),X_im_green(:,:,8),X_im_blue(:,:,8));

figure

ax1 = subplot(1,3,1); imshow(uint8(rgb));

ax2 = subplot(1,3,2); imagesc(X_dif); daspect([1,1,1]);

ax3 = subplot(1,3,3); imshow(Xlabel,color_map);

linkaxes([ax1,ax2,ax3])

%-------------------------------------------

% Cluster Analysis

%-------------------------------------------

%first treshold (a)

%I tried both but I am sure they are both wrong

alfa_i = sqrt( (X_im_blue(:,:,8) - X_im_Med(:,:,3)).^2 + (X_im_red(:,:,8) - X_im_Med(:,:,1)).^2 + (X_im_green(:,:,8) - X_im_Med(:,:,2)).^2);

alfa_i = norm((X_im_blue(:,:,8) - X_im_Med(:,:,3)) + (X_im_red(:,:,8) - X_im_Med(:,:,1)) + (X_im_green(:,:,8) - X_im_Med(:,:,2)));

%I tried to extract only the pixel indexes belonging to a certain cluster but it does not work and then I wouldn't know how to continue

blueband = X_im_blue(:,:,8);

idx = find(Xlabel==1);

pixel_blu_class1 = blueband(idx);