"""Analyze coating layer and similarity measures [#senin]_ [#sim]_.
.. [#senin] Senin, P. (2008).
.. [#sim] https://pypi.org/project/similaritymeasures/
"""
import abc
import dataclasses
from typing import TYPE_CHECKING, Generic, TypeVar
import cv2
import numpy as np
import numpy.typing as npt
from numba import njit # type: ignore
from numba.np.extensions import cross2d # type: ignore
from scipy.optimize import root # type: ignore
from scipy.spatial.distance import cdist # type: ignore
from shapely import LineString, offset_curve # type: ignore
from .cache import attrcache
from .substrate import RectSubstrate, SubstrateBase
if TYPE_CHECKING:
from _typeshed import DataclassInstance
__all__ = [
"CoatingLayerBase",
"CoatingLayer",
"RectLayerShape",
"images_XOR",
"images_ANDXOR",
"equidistant_interpolate",
"parallel_curve",
"acm",
"owp",
"integfrechet_G1",
]
SubstTypeVar = TypeVar("SubstTypeVar", bound=SubstrateBase)
"""Type variable for the substrate type of :class:`CoatingLayerBase`."""
DataTypeVar = TypeVar("DataTypeVar", bound="DataclassInstance")
"""Type variable for :attr:`CoatingLayerBase.DataType`."""
[docs]
class CoatingLayerBase(abc.ABC, Generic[SubstTypeVar, DataTypeVar]):
"""Abstract base class for coating layer object.
Coating layer object stores :class:`SubstrateBase` object and target image, which is
a binary image of coated substrate. The role of coating layer object is to acquire
coating layer region by template matching and analyze its shape.
External API can use the following members to get analysis results of
concrete subclasses.
* :attr:`DataType`: Dataclass type for the analysis result.
* :meth:`analyze`: :attr:`DataType` instance containing analysis result.
* :meth:`draw`: Visualized result.
Arguments:
image: Binary target image.
substrate: Substrate instance storing binary reference image.
tempmatch: Pre-computed template matching result.
External constructor can pass this argument to force the template matching
result. If not passed, :meth:`match_template` performs matching.
"""
[docs]
DataType: type[DataTypeVar]
"""Return type of :attr:`analyze`.
Concrete subclass must assign this attribute with dataclass type.
"""
def __init__(
self,
image: npt.NDArray[np.uint8],
substrate: SubstTypeVar,
*,
tempmatch: tuple[tuple[int, ...], float] | None = None,
):
"""Initialize the instance.
*image* is set to be immutable, and template matching is performed if
*tempmatch* is ``None``.
"""
self._image = image
self._image.setflags(write=False)
self._substrate = substrate
if tempmatch is None:
image = self.image
x0, y0, x1, y1 = self.substrate.reference.templateROI
template = self.substrate.reference.image[y0:y1, x0:x1]
self._tempmatch = self.match_template(image, template)
else:
self._tempmatch = tempmatch
@property
[docs]
def image(self) -> npt.NDArray[np.uint8]:
"""Binary target image.
For immutability, this image is not writable.
"""
return self._image
@property
[docs]
def substrate(self) -> SubstTypeVar:
"""Substrate instance which contains substrate and reference information."""
return self._substrate
@property
[docs]
def tempmatch(self) -> tuple[tuple[int, ...], float]:
"""Template matching location and score."""
return self._tempmatch
[docs]
def match_template(
self, image: npt.NDArray[np.uint8], template: npt.NDArray[np.uint8]
) -> tuple[tuple[int, ...], float]:
"""Perform template matching between *image* and *template*.
Template matching is performed using :func:`cv2.matchTemplate` with
:obj:`cv2.TM_SQDIFF_NORMED`. Subclass may override this method to apply
other algorithm.
Arguments:
image: Binary target image.
template: Binary template image.
"""
res = cv2.matchTemplate(image, template, cv2.TM_SQDIFF_NORMED)
score, _, loc, _ = cv2.minMaxLoc(res)
return (tuple(loc), score)
[docs]
def substrate_point(self) -> npt.NDArray[np.int32]:
"""Upper left point of the substrate image in target image.
Returns:
Coordinates in ``(x, y)``.
"""
temp_point, _ = self.tempmatch
t = self.substrate.reference.templateROI[:2]
s = self.substrate.reference.substrateROI[:2]
return np.array(temp_point, dtype=np.int32) - t + s
@attrcache("_coated_substrate")
[docs]
def coated_substrate(self) -> npt.NDArray[np.bool_]:
"""Coated substrate region.
Returns:
Target image without artifacts, e.g., bath surface.
"""
_, img = cv2.connectedComponents(cv2.bitwise_not(self.image))
x, y = (self.substrate_point() + self.substrate.region_points()).T
return np.isin(img, img[y, x])
@attrcache("_extracted_layer")
@abc.abstractmethod
[docs]
def valid(self) -> bool:
"""Return if the analysis can be performed as expected.
Sometimes, the coating layer instance should be constructed but not analyzed at
all. For example, the coating video may contains frames where the capillary
bridge is not ruptured yet. This method allows analyzer to skip such instances.
"""
@abc.abstractmethod
[docs]
def analyze(self) -> DataTypeVar:
"""Return analysis result as dataclass.
Return type must be :attr:`DataType`.
"""
@abc.abstractmethod
[docs]
def draw(self, *args, **kwargs) -> npt.NDArray[np.uint8]:
"""Return visualization result."""
@dataclasses.dataclass
class CoatingLayerData:
"""Analysis data for :class:`CoatingLayer`."""
pass
[docs]
class CoatingLayer(CoatingLayerBase[SubstrateBase, CoatingLayerData]):
"""Basic implementation of coating layer without any analysis.
Arguments:
image: Binary target image.
substrate: Substrate instance.
tempmatch: Pre-computed template matching result.
Examples:
Construct substrate instance first.
.. plot::
:include-source:
:context: reset
>>> import cv2
>>> from finitedepth import get_sample_path, Reference, Substrate
>>> img = cv2.imread(get_sample_path("ref.png"), cv2.IMREAD_GRAYSCALE)
>>> _, bin = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
>>> ref = Reference(bin, (10, 10, 1250, 200), (100, 100, 1200, 500))
>>> subst = Substrate(ref)
>>> import matplotlib.pyplot as plt #doctest: +SKIP
>>> plt.imshow(subst.draw()) #doctest: +SKIP
Then, construct coating layer instance.
.. plot::
:include-source:
:context: close-figs
>>> from finitedepth import CoatingLayer
>>> img = cv2.imread(get_sample_path("coat.png"), cv2.IMREAD_GRAYSCALE)
>>> _, bin = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
>>> coat = CoatingLayer(bin, subst)
>>> plt.imshow(coat.draw()) #doctest: +SKIP
"""
[docs]
DataType = CoatingLayerData
"""Return :obj:`CoatingLayerData`."""
[docs]
def valid(self) -> bool:
"""Return true, as analysis is not performed at all."""
return True
[docs]
def analyze(self):
"""Return empty :class:`CoatingLayerData`."""
return self.DataType()
[docs]
def draw(
self,
subtraction_mode: str = "none",
layer_color: tuple[int, int, int] = (255, 0, 0),
layer_thickness: int = -1,
) -> npt.NDArray[np.uint8]:
"""Subtract the template match result and paint the coating layer.
Arguments:
subtraction_mode ({`'none', 'template', 'substrate', 'full'`}): Subtraction
mode. `'template'` and `'substrate'` removes overlapping template region
and substrate region, respectively. `'full'` removes both.
layer_color: Layer color for :func:`cv2.drawContours`.
layer_thickness: Layer thickness for :func:`cv2.drawContours`.
"""
image = cv2.cvtColor(self.image, cv2.COLOR_GRAY2RGB)
if subtraction_mode not in ["none", "template", "substrate", "full"]:
raise TypeError("Unrecognized subtraction mode: %s" % subtraction_mode)
if subtraction_mode in ["template", "full"]:
x0, y0, x1, y1 = self.substrate.reference.templateROI
tempImg = self.substrate.reference.image[y0:y1, x0:x1]
h, w = tempImg.shape[:2]
(X0, Y0), _ = self.tempmatch
binImg = self.image[Y0 : Y0 + h, X0 : X0 + w]
mask = images_XOR(~binImg.astype(bool), ~tempImg.astype(bool))
image[Y0 : Y0 + h, X0 : X0 + w][~mask] = 255
if subtraction_mode in ["substrate", "full"]:
x0, y0, x1, y1 = self.substrate.reference.substrateROI
substImg = self.substrate.reference.image[y0:y1, x0:x1]
h, w = substImg.shape[:2]
X0, Y0 = self.substrate_point()
binImg = self.image[Y0 : Y0 + h, X0 : X0 + w]
mask = images_XOR(~binImg.astype(bool), ~substImg.astype(bool))
image[Y0 : Y0 + h, X0 : X0 + w][~mask] = 255
cnts, _ = cv2.findContours(
self.extract_layer().astype(np.uint8),
cv2.RETR_EXTERNAL,
cv2.CHAIN_APPROX_NONE,
)
cv2.drawContours(image, cnts, -1, layer_color, layer_thickness)
return image # type: ignore[return-value]
@dataclasses.dataclass
class RectLayerShapeData:
"""Analysis data for :class:`RectLayerShape`.
Arguments:
LayerLength_Left, LayerLength_Right: Length of the layer on each wall.
Conformality: Conformality of the coating layer.
AverageThickness: Average thickness of the coating layer.
Roughness: Roughness of the coating layer.
MaxThickness_Left, MaxThickness_Bottom, MaxThickness_Right: Regional maximum
thicknesses.
MatchError: Template matching error between ``0`` to ``1``.
``0`` means perfect match.
"""
LayerLength_Left: np.float64
LayerLength_Right: np.float64
Conformality: float
AverageThickness: np.float64
Roughness: float
MaxThickness_Left: np.float64
MaxThickness_Bottom: np.float64
MaxThickness_Right: np.float64
MatchError: float
[docs]
class RectLayerShape(CoatingLayerBase[RectSubstrate, RectLayerShapeData]):
"""Coating layer over rectangular substrate.
Arguments:
image: Binary target image.
substrate: Substrate instance.
opening_ksize: Kernel size for morphological operation.
Elements must be zero or odd number.
reconstruct_radius: Radius of the "safe zone" for noise removal.
Imaginary circles with this radius are drawn on bottom corners of the
substrate. Connected components not passing these circles are regarded as
image artifacts.
roughness_measure: Similarity measure to quantify roughness.
`'DTW'`
Dynamice time warping.
`'SDTW'`
Root mean square of dynamic time warping.
`'IFD'`
Approximated integral Fréchet distance. Requires *ifd_gsize*.
ifd_gsize: Grid size to approximate integral Fréchet distance.
Ignored if *roughness_measure* is not `'IFD'`.
tempmatch: Pre-computed template matching result.
Examples:
Construct substrate instance first.
.. plot::
:include-source:
:context: reset
>>> import cv2
>>> from finitedepth import get_sample_path, Reference, RectSubstrate
>>> img = cv2.imread(get_sample_path("ref.png"), cv2.IMREAD_GRAYSCALE)
>>> _, bin = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
>>> ref = Reference(bin, (10, 10, 1250, 200), (100, 100, 1200, 500))
>>> subst = RectSubstrate(ref, 3.0, 1.0, 0.01)
>>> import matplotlib.pyplot as plt #doctest: +SKIP
>>> plt.imshow(subst.draw()) #doctest: +SKIP
Then, construct coating layer instance.
.. plot::
:include-source:
:context: close-figs
>>> from finitedepth import RectLayerShape
>>> img = cv2.imread(get_sample_path("coat.png"), cv2.IMREAD_GRAYSCALE)
>>> _, bin = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
>>> coat = RectLayerShape(bin, subst, (1, 1), 50, "IFD", 3)
>>> plt.imshow(
... coat.draw(conformality_step=10, roughness_step=5)
... ) #doctest: +SKIP
"""
[docs]
DataType = RectLayerShapeData
"""Return :obj:`RectLayerShapeData`."""
def __init__(
self,
image: npt.NDArray[np.uint8],
substrate: RectSubstrate,
opening_ksize: tuple[int, int],
reconstruct_radius: int,
roughness_measure: str,
ifd_gsize: float | None = None,
*,
tempmatch: tuple[tuple[int, ...], float] | None = None,
):
"""Initialize the instance.
Check whether kernel size are zero or odd and roughness measure is valid.
"""
if not all(i == 0 or (i > 0 and i % 2 == 1) for i in opening_ksize):
raise ValueError("Kernel size must be zero or odd.")
if roughness_measure not in ["DTW", "SDTW", "IFD"]:
raise TypeError(f"Unknown roughness measure: {roughness_measure}")
if roughness_measure == "IFD" and ifd_gsize is None:
raise TypeError("ifd_gsize must be set if roughness_measure is IFD.")
super().__init__(image, substrate, tempmatch=tempmatch)
self._opening_ksize = opening_ksize
self._reconstruct_radius = reconstruct_radius
self._roughness_measure = roughness_measure
self._ifd_gsize = float(ifd_gsize) if ifd_gsize is not None else ifd_gsize
@property
[docs]
def opening_ksize(self) -> tuple[int, int]:
"""Kernel size for morphological operation."""
return self._opening_ksize
@property
[docs]
def reconstruct_radius(self) -> int:
"""Radius of the "safe zone" for noise removal."""
return self._reconstruct_radius
@property
[docs]
def roughness_measure(self) -> str:
"""Similarity measure to quantify roughness."""
return self._roughness_measure
@property
[docs]
def ifd_gsize(self) -> float | None:
"""Grid size to quantify roughness if :attr:`roughness_measure` is IFD."""
return self._ifd_gsize
@attrcache("_layer_contours")
[docs]
def layer_contours(self) -> tuple[npt.NDArray[np.int32], ...]:
"""Find contours of coating layer region.
This method finds external contours of :meth:`~.extract_layer`.
Each contour encloses each discrete region of coating layer.
"""
layer_cnts, _ = cv2.findContours(
self.extract_layer().astype(np.uint8),
cv2.RETR_EXTERNAL,
cv2.CHAIN_APPROX_NONE,
)
return tuple(layer_cnts)
@attrcache("_interfaces")
[docs]
def interfaces(self) -> tuple[npt.NDArray[np.int64], ...]:
"""Find solid-liquid interfaces.
A substrate can have contact with multiple discrete coating layer regions,
and a single coating layer region can have multiple contacts to the substrate.
This method returns indices for :meth:`~finitedepth.PolySubstrateBase.contour`
where solid-liquid interfaces start and stop.
Returns:
tuple of arrays. ``i``-th array represents ``i``-th coating layer region in
:meth:`layer_contours`. Shape of the array is ``(N, 2)``, where ``N`` is the
number of contacts the coating layer region makes. Each column represents
starting and ending indices for the interface interval in substrate contour.
Note:
Each interval describes continuous patch on the substrate contour covered
by the layer. To acquire the interface points, slice :attr:`substrate`'s
:meth:`~finitedepth.PolySubstrateBase.contour` with the indices.
"""
subst_cnt = self.substrate.contour() + self.substrate_point()
ret = []
for layer_cnt in self.layer_contours():
H, W = self.image.shape[:2]
lcnt_img = np.zeros((H, W), dtype=np.uint8)
lcnt_img[layer_cnt[..., 1], layer_cnt[..., 0]] = 255
dilated_lcnt = cv2.dilate(lcnt_img, np.ones((3, 3))).astype(bool)
x, y = subst_cnt.transpose(2, 0, 1)
mask = dilated_lcnt[np.clip(y, 0, H - 1), np.clip(x, 0, W - 1)]
# Find indices of continuous True blocks
idxs = np.where(
np.diff(np.concatenate(([False], mask[:, 0], [False]))) == 1
)[0].reshape(-1, 2)
ret.append(idxs)
return tuple(ret)
@attrcache("_contour")
[docs]
def contour(self) -> npt.NDArray[np.int32]:
"""Contour of the entire coated substrate."""
(cnt,), _ = cv2.findContours(
self.coated_substrate().astype(np.uint8),
cv2.RETR_EXTERNAL,
cv2.CHAIN_APPROX_NONE,
)
return cnt
[docs]
def capbridge_broken(self) -> bool:
"""Check if capillary bridge is ruptured.
As substrate is withdrawn from fluid bath, capillary bridge forms between the
coating layer and bulk fluid and then ruptures.
"""
p0 = self.substrate_point()
_, bl, br, _ = self.substrate.contour()[self.substrate.vertices()]
(B,) = p0 + bl
(C,) = p0 + br
top = np.max([B[1], C[1]])
bot = self.image.shape[0]
if top > bot:
# substrate is located outside of the frame
return False
left = B[0]
right = C[0]
roi_binimg = self.image[top:bot, left:right]
return bool(np.any(np.all(roi_binimg, axis=1)))
[docs]
valid = capbridge_broken
@attrcache("_extracted_layer")
@attrcache("_surface")
[docs]
def surface(self) -> tuple[np.int64, np.int64]:
"""Liquid-gas interface of the coating layer.
Substrate surface exposed to air is considered to be covered by coating layer
with zero thickness.
Returns:
Starting and ending indices for the surface interval in coated substrate
contour.
Note:
To acquire the surface points, slice
:meth:`~finitedepth.PolySubstrateBase.contour` with the indices.
"""
if not self.interfaces():
return (np.int64(-1), np.int64(0))
(i0, i1) = np.sort(np.concatenate(self.interfaces()).flatten())[[0, -1]]
subst_cnt = self.substrate.contour() + self.substrate_point()
endpoints = subst_cnt[[i0, i1]]
vec = self.contour() - endpoints.transpose(1, 0, 2)
(I0, I1) = np.argmin(np.linalg.norm(vec, axis=-1), axis=0)
return (I0, I1 + 1)
@attrcache("_uniform_layer")
@attrcache("_conformality")
@attrcache("_roughness")
[docs]
def roughness(self) -> tuple[float, npt.NDArray[np.float64]]:
"""Similarity-based surface roughness of the coating layer.
Returns:
Roughness between layer surface and uniform layer and its pair of points in
curve space.
"""
I0, I1 = self.surface()
surf = self.contour()[I0:I1]
_, ul = self.uniform_layer()
if surf.size == 0 or ul.size == 0:
return (np.nan, np.empty((0, 2), dtype=np.float64))
ul_len = cv2.arcLength(ul.astype(np.float32), closed=False)
if self.roughness_measure == "DTW":
ul = equidistant_interpolate(ul, int(np.ceil(ul_len)))
dist = cdist(np.squeeze(surf, axis=1), np.squeeze(ul, axis=1))
mat = acm(dist)
path = owp(mat)
roughness = mat[-1, -1] / len(path)
pairs = np.concatenate([surf[path[..., 0]], ul[path[..., 1]]], axis=1)
elif self.roughness_measure == "SDTW":
ul = equidistant_interpolate(ul, int(np.ceil(ul_len)))
dist = cdist(np.squeeze(surf, axis=1), np.squeeze(ul, axis=1))
mat = acm(dist**2)
path = owp(mat)
roughness = np.sqrt(mat[-1, -1] / len(path))
pairs = np.concatenate([surf[path[..., 0]], ul[path[..., 1]]], axis=1)
elif self.roughness_measure == "IFD" and self.ifd_gsize is not None:
surf_len = cv2.arcLength(surf, closed=False)
dist, _, pairs = integfrechet_G1(
np.squeeze(surf, axis=1), np.squeeze(ul, axis=1), self.ifd_gsize
)
roughness = dist / (surf_len + ul_len)
else:
roughness = np.nan
pairs = np.empty((0, 2), dtype=np.float64)
return (float(roughness), pairs)
@attrcache("_max_thickness")
[docs]
def max_thickness(self) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
"""Regional maximum thicknesses.
Coating layer is segmented using :meth:`RectSubstrate.sideline_intersections`.
Points on layer surface and sideline for maximum distance in each region are
found.
Returns:
tuple of two arrays. The first array contains maximum thickness values on
left, bottom, and right region. Value of ``0`` indicates no coating layer on
that region. The second array contains points on layer surface and substrate
lines for the maximum thickness.
Shape of the array is ``(3, 2, 2)``; 1st axis indicates left, bottom
and right region, 2nd axis indicates layer surface and substrate line,
and 3rd axis indicates ``(x, y)`` coordinates.
"""
corners = self.substrate.sideline_intersections() + self.substrate_point()
I0, I1 = self.surface()
surface = self.contour()[I0:I1]
thicknesses, points = [], []
for A, B in zip(corners[:-1], corners[1:]):
AB = B - A
mask = np.cross(AB, surface - A) >= 0
pts = surface[mask]
if pts.size == 0:
thicknesses.append(np.float64(0))
points.append(np.array([[-1, -1], [-1, -1]], np.float64))
else:
Ap = pts - A
Proj = A + AB * (np.dot(Ap, AB) / np.dot(AB, AB))[..., np.newaxis]
dist = np.linalg.norm(Proj - pts, axis=-1)
max_idx = np.argmax(dist)
thicknesses.append(dist[max_idx])
points.append(np.stack([pts[max_idx], Proj[max_idx]]))
return (np.array(thicknesses), np.array(points))
[docs]
def analyze(self):
"""Return :class:`RectLayerShapeData`."""
_, B, C, _ = self.substrate.sideline_intersections() + self.substrate_point()
if not self.interfaces():
LEN_L = LEN_R = np.float64(0)
else:
(i0, i1) = np.sort(np.concatenate(self.interfaces()).flatten())[[0, -1]]
subst_cnt = self.substrate.contour() + self.substrate_point()
pts = subst_cnt[[i0, i1]]
Bp = pts - B
BC = C - B
Proj = B + BC * (np.dot(Bp, BC) / np.dot(BC, BC))[..., np.newaxis]
dists = np.linalg.norm(Proj - pts, axis=-1)
(LEN_L,), (LEN_R,) = dists.astype(np.float64)
C, _ = self.conformality()
AVRGTHCK, _ = self.uniform_layer()
ROUGH, _ = self.roughness()
(THCK_L, THCK_B, THCK_R), _ = self.max_thickness()
_, ERR = self.tempmatch
return self.DataType(
LEN_L,
LEN_R,
C,
AVRGTHCK,
ROUGH,
THCK_L,
THCK_B,
THCK_R,
ERR,
)
[docs]
def draw(
self,
background_mode: str = "image",
subtraction_mode: str = "none",
layer_color: tuple[int, int, int] = (255, 0, 0),
layer_thickness: int = -1,
contactline_color: tuple[int, int, int] = (0, 255, 0),
contactline_thickness: int = 1,
maxthickness_color: tuple[int, int, int] = (0, 255, 0),
maxthickness_thickness: int = 1,
uniformlayer_color: tuple[int, int, int] = (0, 0, 255),
uniformlayer_thickness: int = 1,
conformality_color: tuple[int, int, int] = (0, 0, 255),
conformality_thickness: int = 1,
conformality_step: int = 1,
roughness_color: tuple[int, int, int] = (0, 0, 255),
roughness_thickness: int = 1,
roughness_step: int = 1,
) -> npt.NDArray[np.uint8]:
"""Visualize the analysis result.
#. Draw the substrate with by :class:`PaintMode`.
#. Display the template matching result with :class:`SubtractionMode`.
#. Draw coating layer and contact line.
#. If capillary bridge is broken, draw regional maximum thicknesses,
uniform layer, conformality pairs and roughness pairs.
Arguments:
background_mode ({`'image', 'empty'`}): Determine how background is drawn.
`'image'` draws original background image while `'empty'` draws on
empty frame.
subtraction_mode ({`'none', 'template', 'substrate', 'full'`}): Subtraction
mode. `'template'` and `'substrate'` removes overlapping template region
and substrate region, respectively. `'full'` removes both.
layer_color: Layer contour's color for :func:`cv2.drawContours`.
layer_thickness: Layer contour's thickness for :func:`cv2.drawContours`.
contactline_color: Contact line's color for :func:`cv2.line`.
contactline_thickness: Contact line's thickness for :func:`cv2.line`.
maxthickness_color: Regional maximum thickness line's color for
:func:`cv2.polylines`.
maxthickness_thickness: Regional maximum thickness line's thickness for
:func:`cv2.polylines`.
uniformlayer_color: Imaginary uniform layer's color for
:func:`cv2.polylines`.
uniformlayer_thickness: Imaginary uniform layer's thickness for
:func:`cv2.polylines`.
conformality_color: Conformality pairs' color for :func:`cv2.polylines`.
conformality_thickness: Conformality pairs' thickness for
:func:`cv2.polylines`.
conformality_step: Step size to skip conformality pairs.
roughness_color: Roughness pairs' color for :func:`cv2.polylines`.
roughness_thickness: Roughness pairs' thickness for :func:`cv2.polylines`.
roughness_step: Step size to skip roughness pairs.
"""
if background_mode == "image":
image = self.image
elif background_mode == "empty":
image = np.full(self.image.shape, 255, dtype=np.uint8)
else:
raise TypeError("Unrecognized paint mode: %s" % background_mode)
image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB) # type: ignore[assignment]
if subtraction_mode not in ["none", "template", "substrate", "full"]:
raise TypeError("Unrecognized subtraction mode: %s" % subtraction_mode)
if subtraction_mode in ["template", "full"]:
x0, y0, x1, y1 = self.substrate.reference.templateROI
tempImg = self.substrate.reference.image[y0:y1, x0:x1]
h, w = tempImg.shape[:2]
(X0, Y0), _ = self.tempmatch
binImg = self.image[Y0 : Y0 + h, X0 : X0 + w]
mask = images_XOR(~binImg.astype(bool), ~tempImg.astype(bool))
image[Y0 : Y0 + h, X0 : X0 + w][~mask] = 255
if subtraction_mode in ["substrate", "full"]:
x0, y0, x1, y1 = self.substrate.reference.substrateROI
substImg = self.substrate.reference.image[y0:y1, x0:x1]
h, w = substImg.shape[:2]
X0, Y0 = self.substrate_point()
binImg = self.image[Y0 : Y0 + h, X0 : X0 + w]
mask = images_XOR(~binImg.astype(bool), ~substImg.astype(bool))
image[Y0 : Y0 + h, X0 : X0 + w][~mask] = 255
cv2.drawContours(image, self.layer_contours(), -1, layer_color, layer_thickness)
if len(self.interfaces()) > 0:
(i0, i1) = np.sort(np.concatenate(self.interfaces()).flatten())[[0, -1]]
subst_cnt = self.substrate.contour() + self.substrate_point()
(p0,), (p1,) = subst_cnt[[i0, i1]].astype(np.int32)
cv2.line(image, p0, p1, contactline_color, contactline_thickness)
if not self.capbridge_broken():
return image
lines = []
for dist, pts in zip(*self.max_thickness()):
if dist == 0:
continue
lines.append(pts.astype(np.int32))
cv2.polylines(
image,
lines,
isClosed=False,
color=maxthickness_color,
thickness=maxthickness_thickness,
)
_, points = self.uniform_layer()
cv2.polylines(
image,
[points.astype(np.int32)],
isClosed=False,
color=uniformlayer_color,
thickness=uniformlayer_thickness,
)
if len(self.interfaces()) > 0:
_, pairs = self.conformality()
cv2.polylines(
image,
pairs[::conformality_step],
isClosed=False,
color=conformality_color,
thickness=conformality_thickness,
)
if len(self.interfaces()) > 0:
_, pairs = self.roughness()
cv2.polylines(
image,
pairs.astype(np.int32)[::roughness_step],
isClosed=False,
color=roughness_color,
thickness=roughness_thickness,
)
return image
[docs]
def images_XOR(
img1: npt.NDArray[np.bool_],
img2: npt.NDArray[np.bool_],
point: tuple[int, int] = (0, 0),
) -> npt.NDArray[np.bool_]:
"""Perform subtraction between two images by XOR operation.
Arguments:
img1: Image from which *img2* is subtracted.
img2: Image patch which is subtracted from *img1*.
point: Location in *img1* where *img2* is subtracted.
"""
H, W = img1.shape
h, w = img2.shape
x0, y0 = point
x1, y1 = x0 + w, y0 + h
img1 = img1.copy()
img1_crop = img1[max(y0, 0) : min(y1, H), max(x0, 0) : min(x1, W)]
img2_crop = img2[max(-y0, 0) : min(H - y0, h), max(-x0, 0) : min(W - x0, w)]
img1_crop ^= img2_crop
return img1
[docs]
def images_ANDXOR(
img1: npt.NDArray[np.bool_],
img2: npt.NDArray[np.bool_],
point: tuple[int, int] = (0, 0),
) -> npt.NDArray[np.bool_]:
"""Perform subtraction between two images by AND and XOR operations.
Arguments:
img1: Image from which *img2* is subtracted.
img2: Image patch which is subtracted from *img1*.
point: Location in *img1* where *img2* is subtracted.
"""
H, W = img1.shape
h, w = img2.shape
x0, y0 = point
x1, y1 = x0 + w, y0 + h
img1 = img1.copy()
img1_crop = img1[max(y0, 0) : min(y1, H), max(x0, 0) : min(x1, W)]
img2_crop = img2[max(-y0, 0) : min(H - y0, h), max(-x0, 0) : min(W - x0, w)]
common = img1_crop & img2_crop
img1_crop ^= common
return img1
@njit(cache=True)
[docs]
def equidistant_interpolate(points, n) -> npt.NDArray[np.float64]:
"""Interpolate points with equidistant new points.
Arguments:
points: Points that are interpolated.
The shape must be ``(N, 1, D)`` where ``N`` is the number of points
and ``D`` is the dimension.
n: Number of new points.
Returns:
Interpolated points. If ``N`` is positive number, the shape is ``(n, 1, D)``.
If ``N`` is zero, the shape is ``(n, 0, D)``.
"""
# https://stackoverflow.com/a/19122075
N, _, D = points.shape
if points.size == 0:
return np.empty((n, 0, D), dtype=np.float64)
dist = np.empty((N - 1), dtype=np.float64)
for i in range(N - 1):
dist[i] = np.linalg.norm(points[i + 1, 0, :] - points[i, 0, :])
u = np.empty((N,), dtype=np.float64)
u[0] = 0
u[1:] = np.cumsum(dist)
t = np.linspace(0, u[-1], n)
ret = np.empty((n, 1, D))
for i in range(D):
ret[:, 0, i] = np.interp(t, u, points[:, 0, i])
return ret
[docs]
def parallel_curve(curve: npt.NDArray, dist: float) -> npt.NDArray:
"""Return parallel curve of *curve* with offset distance *dist*.
Arguments:
curve: Vertices of a polyline.
The shape is ``(V, 1, D)``, where ``V`` is the number of vertices and
``D`` is the dimension.
dist: offset distance of the parallel curve.
Returns:
Round-joint parallel curve of shape ``(V, 1, D)``.
"""
if dist == 0:
return curve
ret = offset_curve(LineString(np.squeeze(curve, axis=1)), dist, join_style="round")
return np.array(ret.coords)[:, np.newaxis]
@njit(cache=True)
[docs]
def acm(cm: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
"""Compute accumulated cost matrix from local cost matrix.
Arguments:
cm: Local cost matrix.
Returns:
Accumulated cost matrix. The element at `[-1, -1]` is the total sum along the
optimal path. If *cm* is empty, return value is an empty array.
"""
p, q = cm.shape
ret = np.empty((p, q), dtype=np.float64)
if p == 0 or q == 0:
return ret
ret[0, 0] = cm[0, 0]
for i in range(1, p):
ret[i, 0] = ret[i - 1, 0] + cm[i, 0]
for j in range(1, q):
ret[0, j] = ret[0, j - 1] + cm[0, j]
for i in range(1, p):
for j in range(1, q):
ret[i, j] = min(ret[i - 1, j], ret[i, j - 1], ret[i - 1, j - 1]) + cm[i, j]
return ret
@njit(cache=True)
[docs]
def owp(acm: npt.NDArray[np.float64]) -> npt.NDArray[np.int32]:
"""Compute optimal warping path from accumulated cost matrix.
Arguments:
acm: Accumulated cost matrix.
Returns:
Indices for the two series to get the optimal warping path.
"""
p, q = acm.shape
if p == 0 or q == 0:
return np.empty((0, 2), dtype=np.int32)
path = np.empty((p + q - 1, 2), dtype=np.int32)
path_len = np.int32(0)
i, j = p - 1, q - 1
path[path_len] = [i, j]
path_len += np.int32(1)
while i > 0 or j > 0:
if i == 0:
j -= 1
elif j == 0:
i -= 1
else:
d = min(acm[i - 1, j], acm[i, j - 1], acm[i - 1, j - 1])
if acm[i - 1, j] == d:
i -= 1
elif acm[i, j - 1] == d:
j -= 1
else:
i -= 1
j -= 1
path[path_len] = [i, j]
path_len += np.int32(1)
return path[-(len(path) - path_len + 1) :: -1, :]
[docs]
def integfrechet_G1(
T1: npt.NDArray, T2: npt.NDArray, sigma: float
) -> tuple[float, npt.NDArray[np.int32], npt.NDArray[np.float64]]:
"""Approximate integral Fréchet distance over uniform grid [#mahesh]_.
Arguments:
T1, T2: Polyline vertices. shape=(nodes, 2).
sigma: Grid size.
Returns:
Approximated integral Fréchet distance between two polygonal curves, its
optimal path in parameter space, and its point pairs in curve space.
.. [#mahesh] Maheshwari, A., Sack, J. R., & Scheffer, C. (2018).
"""
seg1 = np.linalg.norm(np.diff(T1, axis=0), axis=1)
seg2 = np.linalg.norm(np.diff(T2, axis=0), axis=1)
return _integfrechet_G1(T1, T2, seg1, seg2, sigma)
@njit(cache=True)
def _integfrechet_G1(T1, T2, seg1, seg2, sigma):
cumsum1 = np.cumsum(seg1)
h = cumsum1[-1]
cumsum2 = np.cumsum(seg2)
b = cumsum2[-1]
i1 = int(np.ceil(h / sigma) + 1)
dy = h / (i1 - 1)
i2 = int(np.ceil(b / sigma) + 1)
dx = b / (i2 - 1)
dists = np.empty((i1, i2, 2), dtype=np.float64) # [..., 0]: +y, [..., 1]: +x
# acm
P, Q = T1[0], T2[0]
dists[0, 0] = [0, 0]
for i in range(1, i1):
dists[i, 0, 0] = dists[i - 1, 0, 0] + _gridedge_integrate(
Q, T1, seg1, cumsum1, dy * (i - 1), dy * i
)
for j in range(1, i2):
dists[0, j, 1] = dists[0, j - 1, 0] + _gridedge_integrate(
P, T2, seg2, cumsum2, dx * (j - 1), dx * j
)
for i in range(1, i1):
p = _polyline_pt(T1, seg1, cumsum1, dy * i)
for j in range(1, i2):
q = _polyline_pt(T2, seg2, cumsum2, dx * j)
dists[i, j, 0] = min(dists[i - 1, j]) + _gridedge_integrate(
q, T1, seg1, cumsum1, dy * (i - 1), dy * i
)
dists[i, j, 1] = min(dists[i, j - 1]) + _gridedge_integrate(
p, T2, seg2, cumsum2, dx * (j - 1), dx * j
)
# owp
path = np.empty((i1 + i2 - 1, 2), dtype=np.int32)
i, j = i1 - 1, i2 - 1
path[i + j] = [i, j]
while i > 0 or j > 0:
if i == 0:
j -= 1
elif j == 0:
i -= 1
else:
d = dists[i, j]
if min(d) == d[0]:
i -= 1
else:
j -= 1
path[i + j] = [i, j]
pairs = np.empty((i1 + i2 - 1, 2, 2), dtype=np.float64)
for i, j in path:
pairs[i + j, 0] = _polyline_pt(T1, seg1, cumsum1, dy * i)
pairs[i + j, 1] = _polyline_pt(T2, seg2, cumsum2, dx * j)
return min(dists[-1, -1]), path, pairs
@njit(cache=True)
def _polyline_pt(poly, seg, cumsum, param):
i = np.searchsorted(cumsum, param)
A, B = poly[i], poly[i + 1]
if i == 0:
t = param / seg[i]
else:
t = (param - cumsum[i - 1]) / seg[i]
return A + (B - A) * t
@njit(cache=True)
def _gridedge_integrate(p, poly, seg, cumsum, param1, param2):
i1, i2 = np.searchsorted(cumsum, param1), np.searchsorted(cumsum, param2)
a = _polyline_pt(poly, seg, cumsum, param1)
b = _polyline_pt(poly, seg, cumsum, param2)
if i1 == i2:
return _dist_integrate(a, b, p)
if i2 > i1:
ret = 0
prev = a.astype(np.float64)
for j in np.arange(i2 - i1):
pt = poly[i1 + j + 1]
ret += _dist_integrate(prev, pt, p)
prev = pt.astype(np.float64)
ret += _dist_integrate(pt, b, p)
return ret
@njit(cache=True)
def _dist_integrate(a, b, p):
ab = np.dot(b - a, b - a)
if ab < 1e-6:
return 0
if np.abs(cross2d(b - a, p - a)) < 1e-6:
t = np.dot(b - a, p - a) / ab
ap = np.dot(p - a, p - a)
bp = np.dot(p - b, p - b)
if t < 0:
return (bp - ap) / 2
elif t > 1:
return (ap - bp) / 2
else:
return (ap + bp) / 2
A = 2 * np.dot(b - a, a - p) / ab
B = np.dot(a - p, a - p) / ab
integ = (
4 * np.sqrt(1 + A + B)
+ 2 * A * (-np.sqrt(B) + np.sqrt(1 + A + B))
- (A**2 - 4 * B)
* np.log((2 + A + 2 * np.sqrt(1 + A + B)) / (A + 2 * np.sqrt(B)))
) / 8
return ab * integ