Histogram matching in Google Earth Engine

matchHistogram.js

/**
 *
 * Histogram matching in Google Earth Engine with S1 SAR images
 * A very small experiment :D
 *
 * Author: Francis Laclé
 * Year: 2019
 *
 */ 

/**
 * Function that returns a normalized image using unitScale()
 * @param image the image object
 * @param band the band to select
 * @param aoi the area of interest for scoping the reducer
 */
function getNormalizedImage(image, band, aoi) {
  var redux = getReducedImage(image, ee.Reducer.minMax(), aoi);
  return image.unitScale(
    ee.Number(redux.get(band + '_min')),
    ee.Number(redux.get(band + '_max'))
  );
}

/**
 * Function wrapper for a reducer
 * @param image the image object
 * @param reducer the reducer to apply, see API docs
 * @param aoi the area of interest for scoping the reducer
 */
function getReducedImage(image, reducer, aoi) {
  return image.reduceRegion({
    reducer: reducer,
    geometry: aoi,
    scale: 10,
    bestEffort: true
  });
}

/**
 * Function to print a histogram in the range of 0..255
 * @param y is the y-value for the chart, in this case the frequencies
 * @param label the label used for the series
 */
function printHistogram(y, label) {
  var x = ee.List.sequence(0, 255);
  var chart = ui.Chart.array.values(y, 0, x)
      .setSeriesNames(label)
      .setOptions({
        title: 'Histogram',
        hAxis: {'title': 'Bin'},
        vAxis: {'title': 'Frequency'},
        pointSize: 3,
        series: {
          0: {color: 'ff0000'}, 
          1: {color: '00ff00'}, 
          2: {color: '0000ff'}
       }
  });
  print(chart);
}

/**
 * Function that returns a histogram array object
 * @param image the image object
 * @param band the band to select
 * @param aoi the area of interest for scoping the reducer
 */
function getHistogram(image, band, aoi) {
  var redux = image.reduceRegion({
    reducer: ee.Reducer.histogram(256, 1),
    geometry: aoi,
    scale: 10,
    maxPixels: 1e9
  });
  var histoDict = ee.Dictionary(redux.get(band)).get('histogram');
  return histoDict;
}

/**
 * Function that returns a cumulative distribution function 
 * as a list object
 * @param image the image object
 * @param band the band to select
 * @param aoi the area of interest for scoping the reducer
 */
function getCDF(image, band, aoi) {
  var histoArr = ee.List(getHistogram(image.multiply(255), band, aoi));
  var len = histoArr.length().subtract(1);
  var nPixels = ee.Array(histoArr).accum(0).toList().get(len);
  var CDF = ee.Array(histoArr).accum(0).divide(nPixels).multiply(len);
  return [CDF.round().toList(), len];
}

/**
 * Function that returns a source image that has been 
 * matched with the target's image histogram using
 * the remap function from Earth Engine
 * @param imageR the source image object
 * @param imageZ the target image object
 * @param bandR the source band to select
 * @param bandZ the target band to select
 * @param aoi the area of interest for scoping the reducer
 */
function matchHistogram(imageR, imageZ, bandR, bandZ, aoi) {
  var cdfR = getCDF(imageR, bandR, aoi);
  var cdfZ = getCDF(imageZ, bandZ, aoi);
  var lenR = cdfR[1];
  var lenZ = cdfZ[1];
  var map = [];
  
  var cdfR0 = cdfR[0].getInfo();
  var cdfZ0 = cdfZ[0].getInfo();
  
  for (var r = 0; r < cdfR0.length; r++) {
    var max = Infinity;
    var min = 0;
    for (var z = 0; z < cdfZ0.length; z++) {
      var diff = Math.pow(cdfZ0[z] - cdfR0[r], 2);
      if (diff < max) {
        max = diff;
        min = z;
      }
    }
    map.push(min);
  }
  
  var seq = ee.List.sequence(0, lenR);
  var imageM = imageR.multiply(lenR).round().remap(seq, map, 0).divide(lenR);
  return imageM.select(['remapped'], [bandR]);
}

// Aruba ~midpoint: -69.96542164802798,12.51783066594331
var poi = ee.Geometry.Point(-69.96542164802798,12.51783066594331);
var bounds = 12000;
var aoi = poi.buffer(bounds).bounds();
var zoom = 12;
Map.centerObject(aoi, zoom);

var img1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterDate('2016-01-01', '2016-06-01')
    .sort('system:time_start', true)
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('resolution', 'H'))
    .filter(ee.Filter.eq('resolution_meters', 10))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
    .select('VV')
    .filterBounds(aoi);
    
var img2 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterDate('2016-06-01', '2016-12-01')
    .sort('system:time_start', true)
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.eq('resolution', 'H'))
    .filter(ee.Filter.eq('resolution_meters', 10))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
    .select('VV')
    .filterBounds(aoi);
    
img1 = getNormalizedImage(img1.mosaic(), 'VV', aoi);
img2 = getNormalizedImage(img2.mosaic(), 'VV', aoi);

var match = matchHistogram(img1, img2, 'VV', 'VV', aoi);

Map.addLayer(img1.clip(aoi), {}, 'img1');
Map.addLayer(img2.clip(aoi), {}, 'img2');
Map.addLayer(match.clip(aoi), {}, 'match');

var vis = {
  min: 0.0,
  max: 1.0,
  bands: ['VV'],
};

print(img1.clip(aoi).getThumbURL(vis));
print(img2.clip(aoi).getThumbURL(vis));
print(match.clip(aoi).getThumbURL(vis));