cloud_mask_algorithms

// Function that selects and returns the correct cloud mask for a given image.
var selectCM = function(func_image, func_CM) {
  var image_date = ee.Date(func_image.get('system:time_start')).format('YYYY-MM-dd').replace('-','','g');
  var cm_band_names = func_CM.bandNames();
  var list = cm_band_names.map(function(instance) {
    return ee.String(instance).rindex(image_date);
  });
  var right_cm = list.indexOf(19);
  var matched_cm_string = ee.String(cm_band_names.get(right_cm));
  var mask = ee.Image(func_CM.select(matched_cm_string)).unmask().remap([0, 1], [1, 0]);
  var func_image_masked = func_image.updateMask(mask);
  return func_image_masked;
}


// Function to simply clip every Image of an Image Collection.
var clipImage = function(clipping_geometry) {
  var wrap = function(image) {
    var clipped_image = image.clip(clipping_geometry);
    return clipped_image;
  };
  return wrap;
};


// Function to simply clip an Image.
var clipSingleImage = function(image, clipping_geometry) {
  return image.clip(clipping_geometry);
};


// Function to compute new CLOUDY_PIXEL_PERCENTAGE-property (CPP) for clipped images.  
var newCPP = function(area_of_interest) {
  var wrap = function(image) {
    var qa60 = image.select(['QA60']);
    qa60 = qa60.unmask();
    // Get all cloudy pixels from QA60-band.
    var clouds = qa60.expression(
      'denseCloud || cirrus',{
        'denseCloud': qa60.eq(1024),
        'cirrus': qa60.eq(2048)
    });
    // Count all cloudy pixels (all pixels with value 1).
    var countCloudyPixels = clouds.reduceRegion({
      reducer: ee.Reducer.sum(),
      geometry: area_of_interest,
      scale : 10
    });
    // Get all cloud-free pixels from QA60-band.
    var noClouds = qa60.expression(
      'cloudless',{
        'cloudless': qa60.eq(0),
    });
    // Count all cloud-free pixels (all pixels with value 1).
    var countNoCloudPixels = noClouds.reduceRegion({
      reducer: ee.Reducer.sum(),
      geometry: area_of_interest,
      scale: 10
    });
    // Compute new ratio (CPP). #cloudy pixels / #all pixels * 100.
    var nrCloudyPixels = ee.Number(countCloudyPixels.get('QA60')).int();
    var nrNoCloudPixels = ee.Number(countNoCloudPixels.get('QA60')).int();
    var allPixels = ee.Number(nrCloudyPixels.add(nrNoCloudPixels));
    var newCPP = ee.Number(nrCloudyPixels.divide(allPixels).multiply(100));
    return image.set({CLOUDY_PIXEL_PERCENTAGE: newCPP});
  };
  return wrap;
};


// Function to eliminate several same images in a row. 
var getRidOfDoubles = function(image, list) {
  var previous = ee.Image(ee.List(list).get(-1));
  var date_0 = ee.Date(previous.get('system:time_start')).format('YYYY-MM-dd');
  var date_1 = ee.Date(ee.Image(image).get('system:time_start')).format('YYYY-MM-dd');
  
  var marked_image = ee.Algorithms.If(ee.Algorithms.IsEqual(date_0, date_1), 
  marked_image = image.set('repetition?', 1), 
  marked_image = image.set('repetition?', 0));
  
  return ee.List(list).add(marked_image);
};


// Function that returns a normalized feature collection together with mean and std. 
var  computeNormalizationFeatureCollection = function(feature_collection, bands_to_normalize) {
  var bands_to_normalize_server = ee.List(bands_to_normalize);
  
  var mean = feature_collection.reduceColumns({
    reducer: ee.Reducer.mean().forEach(bands_to_normalize_server),
    selectors: bands_to_normalize_server
  });
  
  var sd = feature_collection.reduceColumns({
    reducer: ee.Reducer.stdDev().forEach(bands_to_normalize_server),
    selectors: bands_to_normalize_server
  });
  
  var normalizeFeature = function(feature) {
    feature = ee.Feature(feature);
    var values = bands_to_normalize_server.map(function(i) {
      var dummi_number = ee.Number(feature.get(i)).subtract(mean.get(i)).divide(sd.get(i));
      return dummi_number;
    });  
    var dic = ee.Dictionary.fromLists(bands_to_normalize_server, values);
    return feature.setMulti(dic);
  };
  
  var mapped_feature_collection = feature_collection.map(normalizeFeature);
  mapped_feature_collection = mapped_feature_collection.set('mean', mean, 'sd', sd);
  return mapped_feature_collection;
};


// Function that normalizes an Image based on mean and standart deviation values from the training feature.
var applyNormalizationImage = function(img, mean, sd) {
  var normalized_img = ee.ImageCollection.fromImages(
  img.bandNames().map(function(band){
    band = ee.String(band);
    var img_band = img.select(band);
    return img_band.subtract(ee.Number(mean.get(band))).divide(ee.Number(sd.get(band)))
  })).toBands().rename(img.bandNames());
  return normalized_img;
};


// Function to select the clusters. 
var selectClusters = function(img, background_model, result_clustering, n_clusters, bands_thresholds, region_of_interest) {
  
  var n_clusters_list = ee.List.sequence(0, n_clusters.subtract(1));
  
  var bands_norm_difference = bands_thresholds.map(function(bands) {
    var entry = ee.String(bands).cat("_difference");
    return entry;
  });
  
  var img_joined = img.subtract(background_model).select(bands_thresholds, bands_norm_difference).addBands(img.select(bands_thresholds));
  
  var bands_and_difference_bands = bands_thresholds.cat(bands_norm_difference);
  
  var wrap = n_clusters_list.map(function(i) {
    i = ee.Number(i);
    var img_diff_clus = img_joined.updateMask(result_clustering.eq(i)).select(bands_and_difference_bands);
    var clusteri = img_diff_clus.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: region_of_interest,
      bestEffort: true,
      scale: 10,
      tileScale: 4
    });
  
    var clusteri_diff = clusteri.toArray(bands_norm_difference);
    var clusteri_refl = clusteri.toArray(bands_thresholds);
  
    var clusteri_refl_norm = clusteri_refl.multiply(clusteri_refl).reduce(ee.Reducer.mean(), [0]).sqrt().get([0]);
    var clusteri_diff_mean = clusteri_diff.reduce(ee.Reducer.mean(), [0]).get([0]);
    var clusteri_diff_norm = clusteri_diff.multiply(clusteri_diff).reduce(ee.Reducer.mean(), [0]).sqrt().get([0]);

  
    var multitemporal_score_clusteri = ee.Algorithms.If(clusteri_diff_mean.gt(0), clusteri_diff_norm, clusteri_diff_norm.multiply(-1));
    multitemporal_score_clusteri = result_clustering.eq(i).toFloat().multiply(ee.Number(multitemporal_score_clusteri));
    var reflectance_score_clusteri = result_clustering.eq(i).toFloat().multiply(ee.Number(clusteri_refl_norm));
  
    var multitemporal_score = multitemporal_score_clusteri;
    var reflectance_score = reflectance_score_clusteri;
  
    var scores = ee.Dictionary({multitemporalScore: multitemporal_score, reflectanceScore: reflectance_score});
  
    return scores; 
  
  });
  
  return wrap;
};


// Function to merge all Images of an Image Collection into a single Image with 1 band.
var mergeBands = function(image, previous) {
  return ee.Image(previous).add(image);
};


// Function that combines the resulting list from this script (all individual cloud masks) into a single ee.Image with the different 
// list entries beeing the individual layers of that new image.
var cloudMasks2Image = function(cloud_mask, func_image) {
  cloud_mask = ee.Image(cloud_mask);
  func_image = ee.Image(func_image);
  var imageID = ee.String(cloud_mask.get('image_id'));
  func_image = func_image.addBands(cloud_mask.select([0], [ee.String('CM_').cat(imageID)]));
  return func_image;  
};





// Main code for background estimation, invokes the codes above.
var constructCloudMask = function(instance) {
  // Extract the image with clouds, area of interest (convex hulls around homogenous trees), and clipping geometry (study area boundaries)
  // from the dictionary that has been passed to this function. 
  var image =  ee.Image(ee.Dictionary(instance).get('dic_image'));
  image = image.unmask();    
  var area_of_interest = ee.FeatureCollection(ee.Dictionary(instance).get('dic_area_of_interest'));
  var clipping_geometry = ee.FeatureCollection(ee.Dictionary(instance).get('dic_clipping_geometry'));
  // Define important parameters.
  var bands_model = image.bandNames();
  var bands_threshold = ee.List(["B2", "B3", "B4"]);
  var n_clusters = ee.Number(10);
  var threshold_dif_cloud = ee.Number(0.045).multiply(10000);
  var numPixels = ee.Number(1000);
  var threshold_reflectance = ee.Number(0.15).multiply(10000);
  var growing_ratio = ee.Number(2); 
  var imageDate = image.get('system:time_start');
  var imageID = image.id();
  var fullImageID = ee.String('COPERNICUS/S2/').cat(imageID);
  
  // Steps to calculate background model (mean of the three images chronologically before the cloudy image).
  // Extract start and end date for future Image Collection based on cloudy image.
  var end_date = ee.Date(image.get('system:time_start'));
  var start_date_number = ee.Number.parse(ee.Date(image.get('system:time_start')).format('YYYY')).subtract(1);
  var start_date = ee.Date.parse('YYYY', start_date_number);
  var list_dates = ee.List.sequence(1, 2).set(0, start_date).set(1, end_date);
  var image_system_index = image.get('system:index');
  // Construct Image Collection. Start: 1 Year before image aqucisition date. End: acquisition date.
  var image_collection = ee.ImageCollection('COPERNICUS/S2')
  .filterDate(list_dates.get(0), list_dates.get(1))
  .filterBounds(area_of_interest)
  .filter(ee.Filter.neq('system:index', image_system_index))
  .sort('system:time_start');
  //var image_collection = ee.ImageCollection('COPERNICUS/S2')
  //.filterDate(list_dates.get(0), list_dates.get(1))
  //.filterBounds(ee.Geometry.Point([9.426189582664733,48.28392370408561]));
  // Clip constructed image collection to geometry.
  var clipped_image_collection = image_collection.map(clipImage(clipping_geometry));
  // Calculate new CPP value for all images in clipped image collection.
  var clipped_image_collection_newCPP = clipped_image_collection.map(newCPP(clipping_geometry));
  // Filter image collection for CPP less than 10.
  var filtered_clipped_image_collection_newCPP = clipped_image_collection
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 10));
  // Get rid of doubled images. 
  var filtered_clipped_image_collection_newCPP_unique = filtered_clipped_image_collection_newCPP
  .iterate(getRidOfDoubles, ee.List([ee.Image(0).set('system:time_start', 0, 'repetition?', 1)]));
  var filtered_clipped_image_collection_newCPP_unique_1 = ee.ImageCollection(ee.List(filtered_clipped_image_collection_newCPP_unique)
  .filter(ee.Filter.neq('repetition?', 1)));
  // Slice image collection to the last three images. 
  var filtered_clipped_image_collection_newCPP_unique_limit = filtered_clipped_image_collection_newCPP_unique_1
  .limit(3, 'system:time_start', false)
  .sort('system:time_start');
  // Combine pictures using a median filter. Background Model is finished.
  var median_picture = ee.Image(filtered_clipped_image_collection_newCPP_unique_limit
  .median());
  
  // Calculate difference image (cloud image - background model).
  var difference_image = image.subtract(median_picture);

  // K-means clustering.
  var training = difference_image.sample(clipping_geometry, 10, null, null, numPixels, 0, true, 2);
  var normalized_training = computeNormalizationFeatureCollection(training, bands_model);
  var mean = ee.Dictionary(normalized_training.get('mean'));
  var sd = ee.Dictionary(normalized_training.get('sd'));
  normalized_training = normalized_training.set('mean', null, 'sd', null);    // Eventuell wieder rausnehemn... geändert nach dem 20.5
  var clusterer = ee.Clusterer.wekaKMeans(n_clusters).train(normalized_training);
  var difference_image_normalized = applyNormalizationImage(difference_image, mean, sd);
  var result = difference_image_normalized.cluster(clusterer);
  //result = ee.Image(result).set(
  //  'image_id', imageID);
  
  // Calculate multitemporal score and reflectance score.
  var multi_refl_scores = selectClusters(image, median_picture, result, n_clusters, bands_threshold, clipping_geometry);
  var multitemporal_score = ee.List.sequence(0, n_clusters.subtract(1)).map(function(i) {
    return ee.Dictionary(multi_refl_scores.get(i)).get('multitemporalScore');
  });
  var multitemporal_score_IC = ee.ImageCollection(multitemporal_score);
  var multitemporal_score_IMG = ee.Image(multitemporal_score_IC.iterate(mergeBands, ee.Image([0])));
  
  var reflectance_score = ee.List.sequence(0, n_clusters.subtract(1)).map(function(i) {
    return ee.Dictionary(multi_refl_scores.get(i)).get('reflectanceScore');
  });
  var reflectance_score_IC = ee.ImageCollection(reflectance_score);
  var reflectance_score_IMG = ee.Image(reflectance_score_IC.iterate(mergeBands, ee.Image([0])));
 
  // #Apply thresholds.
  var cloud_score_threshold = ee.Image(ee.Algorithms.If(threshold_reflectance.lte(0), 
  cloud_score_threshold = multitemporal_score_IMG.gt(threshold_dif_cloud), 
  cloud_score_threshold = multitemporal_score_IMG.gt(threshold_dif_cloud).multiply(reflectance_score_IMG.gt(threshold_reflectance))));

  // #Apply openings.
  var kernel = ee.Kernel.circle(growing_ratio);
  cloud_score_threshold = cloud_score_threshold.focal_min({kernel: kernel}).focal_max({kernel: kernel});
  
  // #Mask out non-cloud pixels. 
  cloud_score_threshold = cloud_score_threshold.mask(cloud_score_threshold.eq(1));
  
  // #Add metadata, namely imageID.
  cloud_score_threshold = cloud_score_threshold.set(
    'image_id', imageID,
    'full_image_id', fullImageID,
    'sensing_date', imageDate);
  
  return cloud_score_threshold;                                         // Final Cloud Score

};


exports.selectCM = selectCM;
exports.clipImage = clipImage;  
exports.newCPP = newCPP;
exports.getRidOfDoubles = getRidOfDoubles;
exports.computeNormalizationFeatureCollection = computeNormalizationFeatureCollection;
exports.applyNormalizationImage = applyNormalizationImage;
exports.selectClusters = selectClusters;
exports.mergeBands = mergeBands;
exports.cloudMasks2Image = cloudMasks2Image;
exports.constructCloudMask = constructCloudMask;
exports.clipSingleImage = clipSingleImage;