% (c) Rebecca Dell, University of Cambridge, June 2020. 
% Adapted from the work of Williamson et al. (2018) (doi:10.5194/tc-2018-56). 
    % Who calculated depths following Pope et al. (2016) (doi:10.5194/tc-10-15-2016).


function [lake_depths_c] = S2_depth_mask(rinf_red,B4,combined_matrix,filled_combined)

e1 = strel('disk',5); % to dilate lake mask by 5 pixels
e2 = strel('disk',6); % to dilate lake mask by 6 pixels
g_red = 0.7507; % value from Pope and others (2016) 

[combined_mask_labelled,lakenum_c] = bwlabel(combined_matrix,8);
lake_depths_c = zeros(size(rinf_red));

for c = 1:lakenum_c % defines number of times to loop
    
    lmask_c = zeros(size(B4)); % creates lake mask with size as above 
    tmpref_red_c = B4; % tmpref is equal to values of red band 
    % Get the individual lake mask
    lmask_c(combined_mask_labelled == c) = 1;

    % Dilate lake mask by one ring of pixels to find bare-ice albedo 
    lmaskdil_c = imdilate(lmask_c,e1); % e is defined above - structuring element
    lmaskdil_c2 = imdilate(lmask_c,e2); % e is defined above - structuring element
    ledge_c = lmaskdil_c2 - lmaskdil_c; % just get the edge pixels 
    edgepix_red_c = (tmpref_red_c(ledge_c == 1)); % find the red edge pixel refl 
    ad_red_c = nansum(nansum(edgepix_red_c))/sum(sum(~isnan(edgepix_red_c))); % take the mean of these values (ignore NaNs)  

    % mask reflectance arrays 
    tmpref_red_c(lmask_c == 0) = NaN;  
    
    % Calculate Depths
    tmpdepth_red = (log(ad_red_c - rinf_red) - log(tmpref_red_c - rinf_red)) / (g_red);   
    
    tmpdepth_red(tmpdepth_red<0) = NaN;
  
    tmpdepth_c = tmpdepth_red;
    tmpdepth_c(lmask_c == 0) = 0; % if lmask value = 0, then tmpdepth = 0
    
    mean_lake_depth_c = nansum(nansum(tmpdepth_c))/sum(sum(~isnan(tmpdepth_c)));
 
    % add calculated depths to lake depth array
    lake_depths_c = lake_depths_c + tmpdepth_c;
   
    %% Calculate depth of filled holes
 
    clear lmaskdil_c lmaskdil_c2 ledge edgepix_red_c tmpdepth_red tmpref_red_c lmask_c 
    % filled_combined 1's = filled holes
    holes_mask_c = filled_combined;
    
    template_filled_mask_c = zeros(size(filled_combined)); % creates lake mask with size as above
    specific_combined = zeros(size(filled_combined));
    
    % find specific lake 
    specific_combined(combined_mask_labelled == c) = 2;
    % Add the holes mask and specific lake mask
    holes_and_combined = holes_mask_c + specific_combined;
    % Find where there are holes in the specific lake
    template_filled_mask_c(holes_and_combined ==3)=1;
    
    % assign the specific lake holes the mean lake value
    template_filled_mask_c(template_filled_mask_c==1)= mean_lake_depth_c;
    % add the new hole depths to the overall depth mask
    lake_depths_c = lake_depths_c + template_filled_mask_c;
   
   
    clear holes_mask_c template_filled_mask_c specific_combined holes_and_combined template_filled_mask_c 
end

end