Histogram match - Google Earth Engine
/** | |
* # Histogram matching | |
* # Google Earth Engine with S1 SAR | |
* # 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)); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment