filtering module
hydrafloods.filtering
close_binary(img, window=3, neighborhood=None)
Closing morphological filter. Closing is the erosion of the dialiation of values greater than 1.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
binary image to apply closing filter on |
required |
window |
int | ee.Number |
distance in pixels to consider for closing process |
3 |
neighborhood |
int | ee.Number |
size of the neighborhood to perform fast distance trasnform. Smaller values speed up operations however must be greater than window. If no neighborhood is provided, it will be double the window size. default = None |
None |
Returns:
| Type | Description |
|---|---|
ee.Image |
the closed binary image |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def close_binary(img, window=3, neighborhood=None):
"""Closing morphological filter. Closing is the erosion of the dialiation of
values greater than 1.
args:
img (ee.Image): binary image to apply closing filter on
window (int | ee.Number, optional): distance in pixels to consider for closing process
neighborhood (int | ee.Number, optional): size of the neighborhood to perform fast
distance trasnform. Smaller values speed up operations however must be greater
than window. If no neighborhood is provided, it will be double the window size.
default = None
returns:
ee.Image: the closed binary image
"""
if not isinstance(window, ee.Number):
window = ee.Number(window)
if neighborhood is None:
neighborhood = window.multiply(2)
dialation = img.fastDistanceTransform(neighborhood).sqrt().lt(window)
closed = dialation.Not().fastDistanceTransform(neighborhood).sqrt().gt(window)
return closed.updateMask(img.mask())
gamma_map(img, window=7, enl=4.9, keep_bands=['angle'])
Gamma Map speckle filtering algorithm. Algorithm adapted from https://groups.google.com/g/google-earth-engine-developers/c/a9W0Nlrhoq0/m/tnGMC45jAgAJ.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
Earth engine image object. Expects that imagery is a SAR image |
required |
window |
int |
moving window size to apply filter (i.e. a value of 7 == 7x7 window). default = 7 |
7 |
enl |
float |
equivalent number of looks (enl) per pixel from a SAR scan. See https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/resolutions/level-1-ground-range-detected. default = 4.9 |
4.9 |
keep_bands |
list[str] |
list of band names to drop during filtering and include in the result default = ["angle"] |
['angle'] |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered SAR image using the Gamma Map algorithm |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def gamma_map(img, window=7, enl=4.9, keep_bands=["angle"]):
"""Gamma Map speckle filtering algorithm.
Algorithm adapted from https://groups.google.com/g/google-earth-engine-developers/c/a9W0Nlrhoq0/m/tnGMC45jAgAJ.
args:
img (ee.Image): Earth engine image object. Expects that imagery is a SAR image
window (int, optional): moving window size to apply filter (i.e. a value of 7 == 7x7 window). default = 7
enl (float, optional): equivalent number of looks (enl) per pixel from a SAR scan.
See https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/resolutions/level-1-ground-range-detected.
default = 4.9
keep_bands (list[str], optional): list of band names to drop during filtering and include in the result
default = ["angle"]
returns:
ee.Image: filtered SAR image using the Gamma Map algorithm
"""
band_names = img.bandNames()
if keep_bands is not None:
keep_img = img.select(keep_bands)
proc_bands = band_names.removeAll(keep_bands)
else:
proc_bands = band_names
img = img.select(proc_bands)
# Square kernel, window should be odd (typically 3, 5 or 7)
weights = ee.List.repeat(ee.List.repeat(1, window), window)
midPt = (window // 2) + 1 if (window % 2) != 0 else window // 2
# ~~(window/2) does integer division in JavaScript
kernel = ee.Kernel.fixed(window, window, weights, midPt, midPt, False)
# Convert image from dB to natural values
nat_img = geeutils.db_to_power(img)
# Get mean and variance
mean = nat_img.reduceNeighborhood(ee.Reducer.mean(), kernel)
variance = nat_img.reduceNeighborhood(ee.Reducer.variance(), kernel)
# "Pure speckle" threshold
ci = variance.sqrt().divide(mean) # square root of inverse of enl
# If ci <= cu, the kernel lies in a "pure speckle" area -> return simple mean
cu = 1.0 / math.sqrt(enl)
# If cu < ci < cmax the kernel lies in the low textured speckle area -> return the filtered value
cmax = math.sqrt(2.0) * cu
alpha = ee.Image(1.0 + cu * cu).divide(ci.multiply(ci).subtract(cu * cu))
b = alpha.subtract(enl + 1.0)
d = (
mean.multiply(mean)
.multiply(b)
.multiply(b)
.add(alpha.multiply(mean).multiply(nat_img).multiply(4.0 * enl))
)
f = b.multiply(mean).add(d.sqrt()).divide(alpha.multiply(2.0))
caster = ee.Dictionary.fromLists(
proc_bands, ee.List.repeat("float", proc_bands.length())
)
img1 = (
geeutils.power_to_db(mean.updateMask(ci.lte(cu)))
.rename(proc_bands)
.cast(caster)
)
img2 = (
geeutils.power_to_db(f.updateMask(ci.gt(cu)).updateMask(ci.lt(cmax)))
.rename(proc_bands)
.cast(caster)
)
img3 = img.updateMask(ci.gte(cmax)).rename(proc_bands).cast(caster)
# If ci > cmax do not filter at all (i.e. we don't do anything, other then masking)
output = (
ee.ImageCollection([img1, img2, img3])
.reduce(ee.Reducer.firstNonNull())
.rename(proc_bands)
.clip(img.geometry())
)
if keep_bands is not None:
output = output.addBands(keep_img)
# Compose a 3 band image with the mean filtered "pure speckle", the "low textured" filtered and the unfiltered portions
return output
lee_sigma(img, window=9, sigma=0.9, looks=4, tk=7, keep_bands=['angle'])
Lee Sigma speckle filtering algorithm. Implemented from interpreting https://doi.org/10.1109/TGRS.2008.2002881
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
Earth engine image object. Expects that imagery is a SAR image |
required |
window |
int |
moving window size to apply filter (i.e. a value of 9 == 9x9 window). default = 9 |
9 |
sigma |
float |
sigma lookup value from table 1 in paper. default = 0.9 |
0.9 |
looks |
int |
look intensity value from table 1 in paper. default = 4 |
4 |
tk |
int |
threshold value to determine values in window as point targets. default = 7 |
7 |
keep_bands |
list[str] |
list of band names to drop during filtering and include in the result default = ["angle"] |
['angle'] |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered SAR image using the Lee Sigma algorithm |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def lee_sigma(img, window=9, sigma=0.9, looks=4, tk=7, keep_bands=["angle"]):
"""Lee Sigma speckle filtering algorithm.
Implemented from interpreting https://doi.org/10.1109/TGRS.2008.2002881
args:
img (ee.Image): Earth engine image object. Expects that imagery is a SAR image
window (int, optional): moving window size to apply filter (i.e. a value of 9 == 9x9 window). default = 9
sigma (float, optional): sigma lookup value from table 1 in paper. default = 0.9
looks (int, optional): look intensity value from table 1 in paper. default = 4
tk (int, optional): threshold value to determine values in window as point targets. default = 7
keep_bands (list[str], optional): list of band names to drop during filtering and include in the result
default = ["angle"]
returns:
ee.Image: filtered SAR image using the Lee Sigma algorithm
"""
band_names = img.bandNames()
if keep_bands is not None:
keep_img = img.select(keep_bands)
proc_bands = band_names.removeAll(keep_bands)
else:
proc_bands = band_names
img = img.select(proc_bands)
midPt = (window // 2) + 1 if (window % 2) != 0 else window // 2
kernelWeights = ee.List.repeat(ee.List.repeat(1, window), window)
kernel = ee.Kernel.fixed(window, window, kernelWeights, midPt, midPt)
targetWeights = ee.List.repeat(ee.List.repeat(1, 3), 3)
targetkernel = ee.Kernel.fixed(3, 3, targetWeights, 1, 1)
# Lookup table for range and eta values for intensity
sigmaLookup = ee.Dictionary(
{
1: ee.Dictionary(
{
0.5: ee.Dictionary({"A1": 0.436, "A2": 1.92, "η": 0.4057}),
0.6: ee.Dictionary({"A1": 0.343, "A2": 2.21, "η": 0.4954}),
0.7: ee.Dictionary({"A1": 0.254, "A2": 2.582, "η": 0.5911}),
0.8: ee.Dictionary({"A1": 0.168, "A2": 3.094, "η": 0.6966}),
0.9: ee.Dictionary({"A1": 0.084, "A2": 3.941, "η": 0.8191}),
0.95: ee.Dictionary({"A1": 0.043, "A2": 4.840, "η": 0.8599}),
}
),
2: ee.Dictionary(
{
0.5: ee.Dictionary({"A1": 0.582, "A2": 1.584, "η": 0.2763}),
0.6: ee.Dictionary({"A1": 0.501, "A2": 1.755, "η": 0.3388}),
0.7: ee.Dictionary({"A1": 0.418, "A2": 1.972, "η": 0.4062}),
0.8: ee.Dictionary({"A1": 0.327, "A2": 2.260, "η": 0.4819}),
0.9: ee.Dictionary({"A1": 0.221, "A2": 2.744, "η": 0.5699}),
0.95: ee.Dictionary({"A1": 0.152, "A2": 3.206, "η": 0.6254}),
}
),
3: ee.Dictionary(
{
0.5: ee.Dictionary({"A1": 0.652, "A2": 1.458, "η": 0.2222}),
0.6: ee.Dictionary({"A1": 0.580, "A2": 1.586, "η": 0.2736}),
0.7: ee.Dictionary({"A1": 0.505, "A2": 1.751, "η": 0.3280}),
0.8: ee.Dictionary({"A1": 0.419, "A2": 1.865, "η": 0.3892}),
0.9: ee.Dictionary({"A1": 0.313, "A2": 2.320, "η": 0.4624}),
0.95: ee.Dictionary({"A1": 0.238, "A2": 2.656, "η": 0.5084}),
}
),
4: ee.Dictionary(
{
0.5: ee.Dictionary({"A1": 0.694, "A2": 1.385, "η": 0.1921}),
0.6: ee.Dictionary({"A1": 0.630, "A2": 1.495, "η": 0.2348}),
0.7: ee.Dictionary({"A1": 0.560, "A2": 1.627, "η": 0.2825}),
0.8: ee.Dictionary({"A1": 0.480, "A2": 1.804, "η": 0.3354}),
0.9: ee.Dictionary({"A1": 0.378, "A2": 2.094, "η": 0.3991}),
0.95: ee.Dictionary({"A1": 0.302, "A2": 2.360, "η": 0.4391}),
}
),
}
)
# extract data from lookup
looksDict = ee.Dictionary(sigmaLookup.get(ee.String(str(looks))))
sigmaImage = ee.Dictionary(looksDict.get(ee.String(str(sigma)))).toImage()
a1 = sigmaImage.select("A1")
a2 = sigmaImage.select("A2")
aRange = a2.subtract(a1)
eta = sigmaImage.select("η").pow(2)
img = geeutils.db_to_power(img)
# MMSE estimator
mmseMask = img.gte(a1).Or(img.lte(a2))
mmseIn = img.updateMask(mmseMask)
oneImg = ee.Image(1)
z = mmseIn.reduceNeighborhood(ee.Reducer.mean(), kernel, None, True)
varz = mmseIn.reduceNeighborhood(ee.Reducer.variance(), kernel)
varx = (varz.subtract(z.abs().pow(2).multiply(eta))).divide(oneImg.add(eta))
b = varx.divide(varz)
mmse = oneImg.subtract(b).multiply(z.abs()).add(b.multiply(mmseIn))
# workflow
z99 = ee.Dictionary(
img.reduceRegion(
reducer=ee.Reducer.percentile([99], None, 255, 0.001, 1e6),
geometry=img.geometry(),
scale=10,
bestEffort=True,
)
).toImage()
overThresh = img.gte(z99)
K = overThresh.reduceNeighborhood(ee.Reducer.sum(), targetkernel, None, True)
retainPixel = K.gte(tk)
xHat = geeutils.power_to_db(img.updateMask(retainPixel).unmask(mmse))
output = ee.Image(xHat).rename(proc_bands)
if keep_bands is not None:
output = output.addBands(keep_img)
return output
modified_median_zscore(img, fill_img=None)
Outlier detection and filling on complete DEM using the modified z-score and a median filter Method from Iglewicz, B. and Hoaglin, D.C., 1993. How to detect and handle outliers (Vol. 16). Asq Press.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
Earth engine image object to apply filter on |
required |
fill_img |
ee.Image |
Earth engine image object to fill values resulting from filter. If None provided, the fill operation will be on the input image. Default = None |
None |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered image |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def modified_median_zscore(img, fill_img=None):
"""Outlier detection and filling on complete DEM using the modified z-score and a median filter
Method from Iglewicz, B. and Hoaglin, D.C., 1993. How to detect and handle outliers (Vol. 16). Asq Press.
args:
img (ee.Image): Earth engine image object to apply filter on
fill_img (ee.Image): Earth engine image object to fill values resulting from filter.
If None provided, the fill operation will be on the input image. Default = None
returns:
ee.Image: filtered image
"""
kernel = ee.Kernel.fixed(3, 3, [[1, 1, 1], [1, 1, 1], [1, 1, 1]])
kernel_weighted = ee.Kernel.fixed(3, 3, [[1, 1, 1], [1, 0, 1], [1, 1, 1]])
median = img.focal_median(kernel=kernel)
median_weighted = img.focal_median(kernel=kernel_weighted)
diff = img.subtract(median)
mzscore = diff.multiply(0.6745).divide(diff.abs().focal_median(kernel=kernel))
if fill_img:
filled = fill_img.where(mzscore.gt(3.5), median_weighted)
else:
filled = img.where(mzscore.gt(3.5), median_weighted)
return filled
open_binary(img, window=3, neighborhood=None)
Opening morphological filter. Opening is the dilation of the erosion of values greater than 1.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
binary image to apply opening filter on |
required |
window |
int | ee.Number |
distance in pixels to consider for opening process |
3 |
neighborhood |
int | ee.Number |
size of the neighborhood to perform fast distance trasnform. Smaller values speed up operations however must be greater than window. If no nerighborhood is provided, it will be double the window size. default = None |
None |
Returns:
| Type | Description |
|---|---|
ee.Image |
the opened binary image |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def open_binary(img, window=3, neighborhood=None):
"""Opening morphological filter. Opening is the dilation of the erosion of
values greater than 1.
args:
img (ee.Image): binary image to apply opening filter on
window (int | ee.Number, optional): distance in pixels to consider for opening process
neighborhood (int | ee.Number, optional): size of the neighborhood to perform fast
distance trasnform. Smaller values speed up operations however must be greater
than window. If no nerighborhood is provided, it will be double the window size.
default = None
returns:
ee.Image: the opened binary image
"""
if not isinstance(window, ee.Number):
window = ee.Number(window)
if neighborhood is None:
neighborhood = window.multiply(2)
erosion = img.Not().fastDistanceTransform(neighborhood).sqrt().gt(window)
opened = erosion.fastDistanceTransform(neighborhood).sqrt().lt(window)
return opened.updateMask(img.mask())
p_median(img, window=5, keep_bands=['angle'])
P-Median filter for smoothing imagery. Calculates the average from the median along cross and diagnal pixels of a window
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
Earth engine image object to filter |
required |
window |
int |
moving window size to apply filter (i.e. a value of 5 == 5x5 window). default = 5 |
5 |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered image |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def p_median(img, window=5, keep_bands=["angle"]):
"""P-Median filter for smoothing imagery.
Calculates the average from the median along cross and diagnal pixels of a window
args:
img (ee.Image): Earth engine image object to filter
window (int, optional): moving window size to apply filter (i.e. a value of 5 == 5x5 window). default = 5
returns:
ee.Image: filtered image
"""
def _band_filter(bname):
selector = ee.List([bname])
band_img = img.select(selector)
hv_median = band_img.reduceNeighborhood(ee.Reducer.median(), hv_kernel)
diag_median = band_img.reduceNeighborhood(ee.Reducer.median(), diag_kernel)
return ee.Image(ee.Image.cat([hv_median, diag_median]).reduce("mean")).rename(
selector
)
if window % 2 == 0:
window += 1
center_idx = (window - 1) // 2
hv = [
[1 if i == center_idx or j == center_idx else 0 for j in range(window)]
for i in range(window)
]
diag = [
[1 if i == j or i == ((window - 1) - j) else 0 for j in range(window)]
for i in range(window)
]
# method based on ???
band_names = img.bandNames()
hv_weights = ee.List(hv)
diag_weights = ee.List(diag)
hv_kernel = ee.Kernel.fixed(window, window, hv_weights)
diag_kernel = ee.Kernel.fixed(window, window, diag_weights)
reduced_bands = ee.ImageCollection.fromImages(
band_names.map(_band_filter)
).toBands()
return reduced_bands.rename(band_names)
perona_malik(img, n_iters=10, K=3, method=1)
Perona-Malik (anisotropic diffusion) convolution Developed by Gennadii Donchyts see https://groups.google.com/g/google-earth-engine-developers/c/umGlt5qIN1I/m/PD8lsJ7qBAAJ I(n+1, i, j) = I(n, i, j) + lambda * (cN * dN(I) + cS * dS(I) + cE * dE(I), cW * dW(I))
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
img |
ee.Image |
Earth engine image object. Expects that imagery is a SAR image |
required |
n_iters |
int |
Number of interations to apply filter K (int,optional): moving window size to apply filter (i.e. a value of 7 == 7x7 window). default = 3 |
10 |
method |
int |
choose method 1 (default) or 2 |
1 |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered SAR image using the perona malik algorithm |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def perona_malik(img, n_iters=10, K=3, method=1):
"""Perona-Malik (anisotropic diffusion) convolution
Developed by Gennadii Donchyts see https://groups.google.com/g/google-earth-engine-developers/c/umGlt5qIN1I/m/PD8lsJ7qBAAJ
I(n+1, i, j) = I(n, i, j) + lambda * (cN * dN(I) + cS * dS(I) + cE * dE(I), cW * dW(I))
args:
img (ee.Image): Earth engine image object. Expects that imagery is a SAR image
n_iters (int, optional): Number of interations to apply filter
K (int,optional): moving window size to apply filter (i.e. a value of 7 == 7x7 window). default = 3
method (int, optional): choose method 1 (default) or 2
returns:
ee.Image: filtered SAR image using the perona malik algorithm
"""
def _method_1(dI_W, dI_E, dI_N, dI_S):
cW = dI_W.pow(2).multiply(k1).exp()
cE = dI_E.pow(2).multiply(k1).exp()
cN = dI_N.pow(2).multiply(k1).exp()
cS = dI_S.pow(2).multiply(k1).exp()
return cW, cE, cN, cS
def _method_2(dI_W, dI_E, dI_N, dI_S):
cW = one.divide(one.add(dI_W.pow(2).divide(k2)))
cE = one.divide(one.add(dI_E.pow(2).divide(k2)))
cN = one.divide(one.add(dI_N.pow(2).divide(k2)))
cS = one.divide(one.add(dI_S.pow(2).divide(k2)))
return cW, cE, cN, cS
# covnert db to natural units before applying filter
power = geeutils.db_to_power(img)
dxW = ee.Kernel.fixed(3, 3, [[0, 0, 0], [1, -1, 0], [0, 0, 0]])
dxE = ee.Kernel.fixed(3, 3, [[0, 0, 0], [0, -1, 1], [0, 0, 0]])
dyN = ee.Kernel.fixed(3, 3, [[0, 1, 0], [0, -1, 0], [0, 0, 0]])
dyS = ee.Kernel.fixed(3, 3, [[0, 0, 0], [0, -1, 0], [0, 1, 0]])
one = ee.Image.constant(1.0)
l = ee.Image.constant(0.2)
k = ee.Image.constant(K)
k1 = ee.Image.constant(-1.0 / K)
k2 = k.pow(2)
if method == 1:
_method = _method_1
elif method == 2:
_method = _method_2
else:
raise NotImplementedError(
"Could not determine algorithm to apply filter...options for `method` are 1 or 2"
)
for i in range(n_iters):
dI_W = power.convolve(dxW)
dI_E = power.convolve(dxE)
dI_N = power.convolve(dyN)
dI_S = power.convolve(dyS)
cW, cE, cN, cS = _method(dI_W, dI_E, dI_N, dI_S)
power = power.add(
l.multiply(
cN.multiply(dI_N)
.add(cS.multiply(dI_S))
.add(cE.multiply(dI_E))
.add(cW.multiply(dI_W))
)
)
# covnert natural to db units after filter is done
img = geeutils.power_to_db(power)
return img
refined_lee(image, keep_bands=['angle'])
Refined Lee speckle filtering algorithm. Algorithm adapted from https://groups.google.com/g/google-earth-engine-developers/c/ExepnAmP-hQ/m/7e5DnjXXAQAJ
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
image |
ee.Image |
Earth engine image object. Expects that imagery is a SAR image |
required |
Returns:
| Type | Description |
|---|---|
ee.Image |
filtered SAR image using the Refined Lee algorithm |
Source code in hydrafloods/filtering.py
@decorators.keep_names
@decorators.keep_attrs
def refined_lee(image, keep_bands=["angle"]):
"""Refined Lee speckle filtering algorithm.
Algorithm adapted from https://groups.google.com/g/google-earth-engine-developers/c/ExepnAmP-hQ/m/7e5DnjXXAQAJ
args:
image (ee.Image): Earth engine image object. Expects that imagery is a SAR image
returns:
ee.Image: filtered SAR image using the Refined Lee algorithm
"""
# TODO: include keep bands...maybe one-shot filtering if using keep_bands???
def apply_filter(b):
"""Closure function to apply the refined lee algorithm on individual bands"""
b = ee.String(b)
img = power.select(b)
# img must be in natural units, i.e. not in dB!
# Set up 3x3 kernels
weights3 = ee.List.repeat(ee.List.repeat(1, 3), 3)
kernel3 = ee.Kernel.fixed(3, 3, weights3, 1, 1, False)
mean3 = img.reduceNeighborhood(ee.Reducer.mean(), kernel3)
variance3 = img.reduceNeighborhood(ee.Reducer.variance(), kernel3)
# Use a sample of the 3x3 windows inside a 7x7 windows to determine gradients and directions
sample_weights = ee.List(
[
[0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0],
]
)
sample_kernel = ee.Kernel.fixed(7, 7, sample_weights, 3, 3, False)
# Calculate mean and variance for the sampled windows and store as 9 bands
sample_mean = mean3.neighborhoodToBands(sample_kernel)
sample_var = variance3.neighborhoodToBands(sample_kernel)
# Determine the 4 gradients for the sampled windows
gradients = sample_mean.select(1).subtract(sample_mean.select(7)).abs()
gradients = gradients.addBands(
sample_mean.select(6).subtract(sample_mean.select(2)).abs()
)
gradients = gradients.addBands(
sample_mean.select(3).subtract(sample_mean.select(5)).abs()
)
gradients = gradients.addBands(
sample_mean.select(0).subtract(sample_mean.select(8)).abs()
)
# And find the maximum gradient amongst gradient bands
max_gradient = gradients.reduce(ee.Reducer.max())
# Create a mask for band pixels that are the maximum gradient
gradmask = gradients.eq(max_gradient)
# duplicate gradmask bands: each gradient represents 2 directions
gradmask = gradmask.addBands(gradmask)
# Determine the 8 directions
directions = (
sample_mean.select(1)
.subtract(sample_mean.select(4))
.gt(sample_mean.select(4).subtract(sample_mean.select(7)))
.multiply(1)
)
directions = directions.addBands(
sample_mean.select(6)
.subtract(sample_mean.select(4))
.gt(sample_mean.select(4).subtract(sample_mean.select(2)))
.multiply(2)
)
directions = directions.addBands(
sample_mean.select(3)
.subtract(sample_mean.select(4))
.gt(sample_mean.select(4).subtract(sample_mean.select(5)))
.multiply(3)
)
directions = directions.addBands(
sample_mean.select(0)
.subtract(sample_mean.select(4))
.gt(sample_mean.select(4).subtract(sample_mean.select(8)))
.multiply(4)
)
# The next 4 are the not() of the previous 4
directions = directions.addBands(directions.select(0).Not().multiply(5))
directions = directions.addBands(directions.select(1).Not().multiply(6))
directions = directions.addBands(directions.select(2).Not().multiply(7))
directions = directions.addBands(directions.select(3).Not().multiply(8))
# Mask all values that are not 1-8
directions = directions.updateMask(gradmask)
# "collapse" the stack into a singe band image (due to masking, each pixel has just one value (1-8) in it's directional band, and is otherwise masked)
directions = directions.reduce(ee.Reducer.sum())
sample_stats = sample_var.divide(sample_mean.multiply(sample_mean))
# Calculate localNoiseVariance
sigmaV = (
sample_stats.toArray()
.arraySort()
.arraySlice(0, 0, 5)
.arrayReduce(ee.Reducer.mean(), [0])
)
# Set up the 7*7 kernels for directional statistics
rect_weights = ee.List.repeat(ee.List.repeat(0, 7), 3).cat(
ee.List.repeat(ee.List.repeat(1, 7), 4)
)
diag_weights = ee.List(
[
[1, 0, 0, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0, 0],
[1, 1, 1, 0, 0, 0, 0],
[1, 1, 1, 1, 0, 0, 0],
[1, 1, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1],
]
)
rect_kernel = ee.Kernel.fixed(7, 7, rect_weights, 3, 3, False)
diag_kernel = ee.Kernel.fixed(7, 7, diag_weights, 3, 3, False)
# Create stacks for mean and variance using the original kernels. Mask with relevant direction.
dir_mean = img.reduceNeighborhood(ee.Reducer.mean(), rect_kernel).updateMask(
directions.eq(1)
)
dir_var = img.reduceNeighborhood(ee.Reducer.variance(), rect_kernel).updateMask(
directions.eq(1)
)
dir_mean = dir_mean.addBands(
img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(
directions.eq(2)
)
)
dir_var = dir_var.addBands(
img.reduceNeighborhood(ee.Reducer.variance(), diag_kernel).updateMask(
directions.eq(2)
)
)
# and add the bands for rotated kernels
for i in range(1, 4):
dir_mean = dir_mean.addBands(
img.reduceNeighborhood(
ee.Reducer.mean(), rect_kernel.rotate(i)
).updateMask(directions.eq(2 * i + 1))
)
dir_var = dir_var.addBands(
img.reduceNeighborhood(
ee.Reducer.variance(), rect_kernel.rotate(i)
).updateMask(directions.eq(2 * i + 1))
)
dir_mean = dir_mean.addBands(
img.reduceNeighborhood(
ee.Reducer.mean(), diag_kernel.rotate(i)
).updateMask(directions.eq(2 * i + 2))
)
dir_var = dir_var.addBands(
img.reduceNeighborhood(
ee.Reducer.variance(), diag_kernel.rotate(i)
).updateMask(directions.eq(2 * i + 2))
)
# "collapse" the stack into a single band image (due to masking, each pixel has just one value in it's directional band, and is otherwise masked)
dir_mean = dir_mean.reduce(ee.Reducer.sum())
dir_var = dir_var.reduce(ee.Reducer.sum())
# A finally generate the filtered value
varX = dir_var.subtract(dir_mean.multiply(dir_mean).multiply(sigmaV)).divide(
sigmaV.add(1.0)
)
b = varX.divide(dir_var)
# return multi-band image band from array
return (
dir_mean.add(b.multiply(img.subtract(dir_mean)))
.arrayProject([0])
.arrayFlatten([["sum"]])
.float()
)
band_names = image.bandNames()
if keep_bands is not None:
keep_img = image.select(keep_bands)
proc_bands = band_names.removeAll(keep_bands)
else:
proc_bands = band_names
image = image.select(proc_bands)
power = geeutils.db_to_power(image)
result = ee.ImageCollection(proc_bands.map(apply_filter)).toBands()
output = geeutils.power_to_db(ee.Image(result)).rename(proc_bands)
if keep_bands is not None:
output = output.addBands(keep_img)
return output