Otsu's method combined with water index to calculate the time-series data?

Question

This question is a continuation of Calculate water frequency in Google Earth Engine based on time series?

The previous question used a set of criteria (the numerical relationship like "larger than", "less than" etc between some indices to extract the water pixels. Then applied a decision tree to get the time series map of water frequency: the map shows the frequency of each pixel appear as water among a collection Landsat images of the whole year.

To get the better result, I slightly revised the code and used a new index: AWEIs, to extract the water pixels.

The difference is, in the previous question, the pixels are filtered by numerical relationship, so the output is simply "Yes" or "No", so the binary map is easy to generated.

However, now I use an index, which is a formula and the output are numbers. So, if I want to get a water-land binary map, I need to find a cutoff value to set the threshold.

To easy demonstrate it, here I attached my revised code, which uses the index of AWEIs:

//Choose Area
var region = ee.Geometry.Polygon({
  coords: [[[-80.79, 24.98], [-80.08, 26.43], [-80.71, 26.16],[-81.07, 26.20]]],
  geodesic: false
});

//Add image collection
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;

var non_mosaic = landsat8 .filterBounds(region) .filterDate('2017-01-01', '2017-12-31');

var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out cloudy shadow pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
  // Return an image masking out cloudy areas.
};

// A function to mask out cloud pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 5,5, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
};

//Will be used for frequency calculaton
var clouds_free = non_mosaic.map(maskClouds);

//display the mosaicked origional image and cloud free image
var mosaic_free = non_mosaic.map(maskClouds).median();
var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};
Map.addLayer(non_mosaic, visParams, 'With clouds'); 
Map.addLayer(mosaic_free, visParams, 'Cloud free');

//design the decision tree for water and vegetation coverage detection
var decision_tree = function(image){
  var mndwi = image.normalizedDifference(['B3', 'B6']);
  var evi = image.expression(
    '(NIR - RED)/(NIR+6*RED-7.5*BLUE+1)*2.5', {
      'BLUE': image.select('B2'),
      'RED': image.select('B4'),
      'NIR': image.select('B5'),
    });
  var aweis = image.expression(
    'BLUE+2.5*GREEN-1.5*(NIR+SWIR1)-0.25*SWIR2', {
      'BLUE' : image.select('B2'),
      'GREEN': image.select('B3'),
      'SWIR1': image.select('B6'),
      'NIR': image.select('B5'),
      'SWIR2': image.select('B7'),
});
var msavi = mosaic_free.expression(
    '(2*NIR+1-((2*NIR+1)**2-8*(NIR-RED))**(1/2))/2', {
      'RED' : mosaic_free.select('B4'),
      'NIR': mosaic_free.select('B5'),
});
  var ndvi = image.normalizedDifference(['B5', 'B4']);
  var lswi = image.normalizedDifference(['B5', 'B6']);
  var water = aweis.gte(0).rename('water');
  var vegetation = evi.gte(0.1).and (ndvi.gte(0.2)).and (lswi.gt(0)).rename('vegetation');
  return water.addBands(vegetation).copyProperties(image);  
};

var result = clouds_free.map(decision_tree);
//calculate the frequency
var freq = result.sum().divide(result.count());

var VisParamWar = {"bands":["water"],"min":0,"max":1, palette: ['0000ff', '00ffff','ffffff', 'ffff00', 'ff0000']};
var VisParamVeg = {"bands":["vegetation"],"min":0,"max":1, palette: ['0000ff', '00ffff', 'ffffff', 'ffff00', 'ff0000']};

//display the frequencies
Map.setCenter(-80.702898,25.679568, 8);
Map.addLayer(freq,VisParamWar,'Water Frequency');
Map.addLayer(freq,VisParamVeg,'Vegetation Frequency');

The above code works well. As you can see, I selected 0 as the threshold of AWEIs. The pixels with the values equal or larger than 0 are considered as water pixels.

Now, instead of the constant value, I want to use Otsu's method to find the optimized threshold and apply it to the water frequency calculation. Here I share my code as follow: