Skip to content

Commit

Permalink
Switch morphology tasks to use cv2 instead of skimage (#732)
Browse files Browse the repository at this point in the history 
* update tests for erosion

* add morphological filter regression tasks

* remove usage of skimage in morphology

* update tests

---------

Co-authored-by: Ziga Luksic <[email protected]>
  • Loading branch information
@zigaLuksic
Matic Lubej and zigaLuksic authored Aug 30, 2023
1 parent 98cecf1 commit 6037f4d
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 68 deletions.
4 changes: 2 additions & 2 deletions eolearn/core/utils/testing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@

from ..constants import FeatureType
from ..eodata import EOPatch
from ..types import Feature
from ..types import FeaturesSpecification
from ..utils.parsing import FeatureParser

DEFAULT_BBOX = BBox((0, 0, 100, 100), crs=CRS("EPSG:32633"))
Expand All @@ -46,7 +46,7 @@ def __post_init__(self) -> None:


def generate_eopatch(
features: list[Feature] | None = None,
features: FeaturesSpecification | None = None,
bbox: BBox = DEFAULT_BBOX,
timestamps: list[dt.datetime] | None = None,
seed: int = 42,
Expand Down
47 changes: 17 additions & 30 deletions eolearn/geometry/morphology.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,11 +10,10 @@

import itertools as it
from enum import Enum
from typing import Callable
from functools import partial

import cv2
import numpy as np
import skimage.filters.rank
import skimage.morphology

from eolearn.core import EOPatch, EOTask, MapFeatureTask
from eolearn.core.types import FeaturesSpecification, SingleFeatureSpec
Expand Down Expand Up @@ -44,7 +43,7 @@ def __init__(
parsed_mask_feature = parse_renamed_feature(mask_feature, allowed_feature_types=lambda fty: fty.is_array())

self.mask_type, self.mask_name, self.new_mask_name = parsed_mask_feature
self.disk = skimage.morphology.disk(disk_radius)
self.disk = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (disk_radius, disk_radius))
self.erode_labels = erode_labels
self.no_data_label = no_data_label

Expand All @@ -57,7 +56,7 @@ def execute(self, eopatch: EOPatch) -> EOPatch:
erode_labels = set(erode_labels) - {self.no_data_label}
other_labels = set(all_labels) - set(erode_labels) - {self.no_data_label}

eroded_masks = [skimage.morphology.binary_erosion(feature_array == label, self.disk) for label in erode_labels]
eroded_masks = [cv2.erode((feature_array == label).astype(np.uint8), self.disk) for label in erode_labels]
other_masks = [feature_array == label for label in other_labels]

merged_mask = np.logical_or.reduce(eroded_masks + other_masks, axis=0)
Expand All @@ -75,20 +74,18 @@ class MorphologicalOperations(Enum):
CLOSING = "closing"
DILATION = "dilation"
EROSION = "erosion"
MEDIAN = "median"

@classmethod
def get_operation(cls, morph_type: MorphologicalOperations) -> Callable:
def get_operation(cls, morph_type: MorphologicalOperations) -> int:
"""Maps morphological operation type to function
:param morph_type: Morphological operation type
"""
return {
cls.OPENING: skimage.morphology.opening,
cls.CLOSING: skimage.morphology.closing,
cls.DILATION: skimage.morphology.dilation,
cls.EROSION: skimage.morphology.erosion,
cls.MEDIAN: skimage.filters.rank.median,
cls.OPENING: cv2.MORPH_OPEN,
cls.CLOSING: cv2.MORPH_CLOSE,
cls.DILATION: cv2.MORPH_DILATE,
cls.EROSION: cv2.MORPH_ERODE,
}[morph_type]


Expand All @@ -103,15 +100,7 @@ def get_disk(radius: int) -> np.ndarray:
:param radius: Radius of disk
:return: The structuring element where elements of the neighborhood are 1 and 0 otherwise.
"""
return skimage.morphology.disk(radius)

@staticmethod
def get_diamond(radius: int) -> np.ndarray:
"""
:param radius: Radius of diamond
:return: The structuring element where elements of the neighborhood are 1 and 0 otherwise.
"""
return skimage.morphology.diamond(radius)
return cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (radius, radius))

@staticmethod
def get_rectangle(width: int, height: int) -> np.ndarray:
Expand All @@ -120,15 +109,15 @@ def get_rectangle(width: int, height: int) -> np.ndarray:
:param height: Height of rectangle
:return: A structuring element consisting only of ones, i.e. every pixel belongs to the neighborhood.
"""
return skimage.morphology.rectangle(width, height)
return cv2.getStructuringElement(cv2.MORPH_RECT, (height, width))

@staticmethod
def get_square(width: int) -> np.ndarray:
"""
:param width: Size of square
:return: A structuring element consisting only of ones, i.e. every pixel belongs to the neighborhood.
"""
return skimage.morphology.square(width)
return cv2.getStructuringElement(cv2.MORPH_RECT, (width, width))


class MorphologicalFilterTask(MapFeatureTask):
Expand All @@ -139,7 +128,7 @@ def __init__(
input_features: FeaturesSpecification,
output_features: FeaturesSpecification | None = None,
*,
morph_operation: MorphologicalOperations | Callable,
morph_operation: MorphologicalOperations,
struct_elem: np.ndarray | None = None,
):
"""
Expand All @@ -153,21 +142,19 @@ def __init__(
output_features = input_features
super().__init__(input_features, output_features)

if isinstance(morph_operation, MorphologicalOperations):
self.morph_operation = MorphologicalOperations.get_operation(morph_operation)
else:
self.morph_operation = morph_operation
self.morph_operation = MorphologicalOperations.get_operation(morph_operation)
self.struct_elem = struct_elem

def map_method(self, feature: np.ndarray) -> np.ndarray:
"""Applies the morphological operation to a raster feature."""
feature = feature.copy()
morph_func = partial(cv2.morphologyEx, kernel=self.struct_elem, op=self.morph_operation)
if feature.ndim == 3:
for channel in range(feature.shape[2]):
feature[..., channel] = self.morph_operation(feature[..., channel], self.struct_elem)
feature[..., channel] = morph_func(feature[..., channel])
elif feature.ndim == 4:
for time, channel in it.product(range(feature.shape[0]), range(feature.shape[3])):
feature[time, ..., channel] = self.morph_operation(feature[time, ..., channel], self.struct_elem)
feature[time, ..., channel] = morph_func(feature[time, ..., channel])
else:
raise ValueError(f"Invalid number of dimensions: {feature.ndim}")

Expand Down
82 changes: 46 additions & 36 deletions tests/geometry/test_morphology.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,8 @@
import pytest
from numpy.testing import assert_array_equal

from sentinelhub import CRS, BBox

from eolearn.core import EOPatch, FeatureType
from eolearn.core.utils.testing import PatchGeneratorConfig, generate_eopatch
from eolearn.geometry import ErosionTask, MorphologicalFilterTask, MorphologicalOperations, MorphologicalStructFactory

CLASSES = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
Expand All @@ -21,65 +20,76 @@
# ruff: noqa: NPY002


@pytest.fixture(name="patch")
def patch_fixture() -> EOPatch:
config = PatchGeneratorConfig(max_integer_value=10, raster_shape=(50, 100), depth_range=(3, 4))
patch = generate_eopatch([MASK_FEATURE, MASK_TIMELESS_FEATURE], config=config)
for feat in [MASK_FEATURE, MASK_TIMELESS_FEATURE]:
patch[feat] = patch[feat].astype(np.uint8)
return patch


@pytest.mark.parametrize("invalid_input", [None, 0, "a"])
def test_erosion_value_error(invalid_input):
with pytest.raises(ValueError):
ErosionTask((FeatureType.MASK_TIMELESS, "LULC", "TEST"), disk_radius=invalid_input)


def test_erosion_full(test_eopatch):
mask_before = test_eopatch.mask_timeless["LULC"].copy()

erosion_task = ErosionTask((FeatureType.MASK_TIMELESS, "LULC", "LULC_ERODED"), 1)
erosion_task = ErosionTask((FeatureType.MASK_TIMELESS, "LULC", "LULC_ERODED"), disk_radius=3)
eopatch = erosion_task.execute(test_eopatch)

mask_after = eopatch.mask_timeless["LULC_ERODED"].copy()
mask_after = eopatch.mask_timeless["LULC_ERODED"]

assert not np.all(mask_before == mask_after)

for label in CLASSES:
if label == 0:
assert np.sum(mask_after == label) >= np.sum(mask_before == label), "Error in the erosion process"
else:
assert np.sum(mask_after == label) <= np.sum(mask_before == label), "Error in the erosion process"
assert_array_equal(np.unique(mask_after, return_counts=True)[1], [1942, 6950, 1069, 87, 52])


def test_erosion_partial(test_eopatch):
mask_before = test_eopatch.mask_timeless["LULC"].copy()

# skip forest and artificial surface
specific_labels = [0, 1, 3, 4]
erosion_task = ErosionTask(
mask_feature=(FeatureType.MASK_TIMELESS, "LULC", "LULC_ERODED"), disk_radius=1, erode_labels=specific_labels
mask_feature=(FeatureType.MASK_TIMELESS, "LULC", "LULC_ERODED"), disk_radius=3, erode_labels=specific_labels
)
eopatch = erosion_task.execute(test_eopatch)

mask_after = eopatch.mask_timeless["LULC_ERODED"].copy()
mask_after = eopatch.mask_timeless["LULC_ERODED"]

assert not np.all(mask_before == mask_after)
assert_array_equal(np.unique(mask_after, return_counts=True)[1], [1145, 7601, 1069, 87, 198])

for label in CLASSES:
if label == 0:
assert np.sum(mask_after == label) >= np.sum(mask_before == label), "Error in the erosion process"
elif label in specific_labels:
assert np.sum(mask_after == label) <= np.sum(mask_before == label), "Error in the erosion process"
else:
assert_array_equal(mask_after == label, mask_before == label, err_msg="Error in the erosion process")


@pytest.mark.parametrize("morph_operation", MorphologicalOperations)
@pytest.mark.parametrize(
"struct_element", [None, MorphologicalStructFactory.get_disk(5), MorphologicalStructFactory.get_rectangle(5, 6)]
("morph_operation", "struct_element", "mask_counts", "mask_timeless_counts"),
[
(
MorphologicalOperations.DILATION,
None,
[6, 34, 172, 768, 2491, 7405, 19212, 44912],
[1, 2, 16, 104, 466, 1490, 3870, 9051],
),
(MorphologicalOperations.EROSION, MorphologicalStructFactory.get_disk(11), [74957, 42, 1], [14994, 6]),
(MorphologicalOperations.OPENING, MorphologicalStructFactory.get_disk(11), [73899, 1051, 50], [14837, 163]),
(MorphologicalOperations.CLOSING, MorphologicalStructFactory.get_disk(11), [770, 74230], [425, 14575]),
(
MorphologicalOperations.OPENING,
MorphologicalStructFactory.get_rectangle(5, 6),
[48468, 24223, 2125, 169, 15],
[10146, 4425, 417, 3, 9],
),
(
MorphologicalOperations.DILATION,
MorphologicalStructFactory.get_rectangle(5, 6),
[2, 19, 198, 3929, 70852],
[32, 743, 14225],
),
],
)
def test_morphological_filter(morph_operation, struct_element):
eopatch = EOPatch(bbox=BBox((0, 0, 1, 1), CRS(3857)), timestamps=["2015-7-7"] * 10)
eopatch[MASK_FEATURE] = np.random.randint(20, size=(10, 100, 100, 3), dtype=np.uint8)
eopatch[MASK_TIMELESS_FEATURE] = np.random.randint(20, 50, size=(100, 100, 5), dtype=np.uint8)

def test_morphological_filter(patch, morph_operation, struct_element, mask_counts, mask_timeless_counts):
task = MorphologicalFilterTask(
[MASK_FEATURE, MASK_TIMELESS_FEATURE], morph_operation=morph_operation, struct_elem=struct_element
)
task.execute(eopatch)
task.execute(patch)

assert eopatch[MASK_FEATURE].shape == (10, 100, 100, 3)
assert eopatch[MASK_TIMELESS_FEATURE].shape == (100, 100, 5)
assert patch[MASK_FEATURE].shape == (5, 50, 100, 3)
assert patch[MASK_TIMELESS_FEATURE].shape == (50, 100, 3)
assert_array_equal(np.unique(patch[MASK_FEATURE], return_counts=True)[1], mask_counts)
assert_array_equal(np.unique(patch[MASK_TIMELESS_FEATURE], return_counts=True)[1], mask_timeless_counts)

0 comments on commit 6037f4d

Please sign in to comment.