"""Finite depth dip coating uniformity measurement using integral Fréchet distance.
To reproduce the examples, run the following code first::
import cv2
from finitedepth import *
from finitedepth_ifd import *
"""
import abc
import dataclasses
from functools import partial
import cv2
import numpy as np
from curvesimilarities import ifd_owp # type: ignore[import-untyped]
from curvesimilarities.util import sample_polyline # type: ignore[import-untyped]
from finitedepth import CoatingLayerBase
from finitedepth.cache import attrcache
from finitedepth.coatinglayer import parallel_curve
from scipy.optimize import root # type: ignore
__all__ = [
"IfdRoughnessBase",
"RectIfdRoughness",
"RectIfdRoughnessData",
]
[docs]
class IfdRoughnessBase(CoatingLayerBase):
"""Base class to measure the coating layer roughness with integral Fréchet distance.
The :class:`IfdRoughnessBase` generalizes the :math:`R_q` roughness [#]_ into
arbitrary geometries by computing the quadratic mean of the integral Fréchet
distance. See :meth:`roughness` for more information.
Parameters
----------
image, substrate
See :class:`CoatingLayerBase <finitedepth.CoatingLayerBase>`.
delta : double
The maximum distance between the Steiner points to compute the roughness.
Refer to :meth:`roughness` for more explanation.
Other Parameters
----------------
tempmatch : tuple, optional
See :class:`CoatingLayerBase <finitedepth.CoatingLayerBase>`.
References
----------
.. [#] https://en.wikipedia.org/wiki/Surface_roughness
"""
def __init__(self, image, substrate, delta, *, tempmatch=None):
if not isinstance(delta, float):
raise TypeError("delta must be a double-precision float.")
if not delta > 0:
raise TypeError("delta must be a positive number.")
super().__init__(image, substrate, tempmatch=tempmatch)
self.delta = delta
[docs]
@abc.abstractmethod
def surface(self):
"""Coating layer surface points.
Returns
-------
ndarray
An :math:`N` by :math:`2` array containing the :math:`xy`-coordinates
of :math:`N` points which constitute the coating layer surface profile.
"""
...
[docs]
@attrcache("_roughness")
def roughness(self):
"""Surface roughness of the coating layer.
Roughness is similarity between :meth:`surface` and :meth:`uniform_layer`.
Here, we choose quadratic average Fréchet distance as the similarity.
The :attr:`delta` attribute determines the approximation accuracy. Refer to the
See Also section for more details.
Returns
-------
roughness : double
Roughness value.
path : ndarray
An :math:`P` by :math:`2` array representing the optimal warping path
in the parameter space.
See Also
--------
curvesimilarities.averagefrechet.qafd : Quadratic average Fréchet distance.
"""
squared_ifd, path = ifd_owp(
self.surface(), self.uniform_layer(), self.delta, dist="squared_euclidean"
)
roughness = np.sqrt(squared_ifd / np.sum(path[-1]))
return float(roughness), path
[docs]
@dataclasses.dataclass
class RectIfdRoughnessData:
"""Analysis data for :class:`RectIfdRoughness`.
Attributes
----------
AverageThickness : double
Average thickness of the coating layer.
Roughness : double
Coating layer roughness.
"""
AverageThickness: float
Roughness: float
[docs]
class RectIfdRoughness(IfdRoughnessBase):
"""Measure coating layer surface roughness over rectangular substrate.
Parameters
----------
image
See :class:`CoatingLayerBase <finitedepth.CoatingLayerBase>`.
substrate : :class:`RectSubstrate <finitedepth.RectSubstrate>`.
Substrate instance.
opening_ksize : tuple of int
Kernel size for morphological opening operation. Must be zero or odd.
reconstruct_radius : int
Radius of the safe zone for noise removal.
Two imaginary circles with this radius are drawn on bottom corners of the
substrate. When extracting the coating layer, connected components not spanning
over any of these circles are regarded as noise.
Other Parameters
----------------
tempmatch : tuple, optional
See :class:`CoatingLayerBase <finitedepth.CoatingLayerBase>`.
Examples
--------
.. note::
For every example in this class, the following code is assumed to be run before.
Construct the substrate instance first.
>>> ref_img = cv2.imread(get_sample_path("ref.png"), cv2.IMREAD_GRAYSCALE)
>>> ref = Reference(
... cv2.threshold(ref_img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1],
... (10, 10, 1250, 200),
... (100, 100, 1200, 500),
... )
>>> subst = RectSubstrate(ref, 3.0, 1.0, 0.01)
Construct the coating layer instance.
>>> target_img = cv2.imread(get_sample_path("coat.png"), cv2.IMREAD_GRAYSCALE)
>>> coat = RectIfdRoughness(
... cv2.threshold(target_img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1],
... subst,
... 5.0,
... (1, 1),
... 50,
... )
Analyze and visualize the coating layer.
>>> coat.analyze()
RectIfdRoughnessData(AverageThickness=50.25..., Roughness=44.91...)
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> plt.imshow(coat.draw()) # doctest: +SKIP
"""
DataType = RectIfdRoughnessData
def __init__(
self,
image,
substrate,
delta,
opening_ksize,
reconstruct_radius,
*,
tempmatch=None,
):
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 reconstruct_radius < 0:
raise ValueError("Reconstruct radius must be zero or positive.")
super().__init__(image, substrate, delta, tempmatch=tempmatch)
self.opening_ksize = opening_ksize
self.reconstruct_radius = reconstruct_radius
[docs]
def valid(self):
"""Check if the coating layer is valid.
The coating layer is invalid if the capillary bridge is not ruptured.
Returns
-------
bool
"""
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]
def substrate_contour(self):
"""Return :attr:`substrate`'s contour in :attr:`image`.
Returns
-------
ndarray
Array of substrate contour points.
"""
return self.substrate.contour() + self.substrate_point()
[docs]
@attrcache("_interface_indices")
def interface_indices(self):
"""Return indices of the substrate contour for the solid-liquid interface.
The interface points can be retrieved by slicing the substrate contour with
the indices.
Returns
-------
ndarray
Starting and ending indices for the solid-liquid interface, empty if the
interface does not exist.
See Also
--------
substrate_contour : The substrate contour which can be sliced.
Examples
--------
.. only:: doctest
>>> coat = getfixture('RectIfdRoughness_setup')
>>> i0, i1 = coat.interface_indices()
>>> interface = coat.substrate_contour()[i0:i1]
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> plt.imshow(coat.image, cmap="gray") # doctest: +SKIP
>>> plt.plot(*interface.transpose(2, 0, 1)) # doctest: +SKIP
Notes
-----
The interface is detected by finding the points on the substrate contour which
are adjacent to the points in :meth:`extract_layer`.
"""
layer_dil = cv2.dilate(self.extract_layer().astype(np.uint8), np.ones((3, 3)))
x, y = self.substrate_contour().transpose(2, 0, 1)
H, W = self.image.shape[:2]
mask = layer_dil[np.clip(y, 0, H - 1), np.clip(x, 0, W - 1)]
idx = np.nonzero(mask[:, 0])[0]
if len(idx) > 0:
idx = idx[[0, -1]]
return idx
[docs]
@attrcache("_surface")
def surface(self):
"""See :meth:`IfdRoughnessBase.surface`.
Examples
--------
>>> import matplotlib.pyplot as plt # doctest: +SKIP
>>> plt.imshow(coat.image, cmap="gray") # doctest: +SKIP
>>> plt.plot(*coat.surface().T, color="tab:red") # doctest: +SKIP
"""
idxs = self.interface_indices()
if len(idxs) == 0:
return np.empty((0, 2), dtype=np.int32)
boundary_pts = self.substrate_contour()[idxs]
(cnt,), _ = cv2.findContours(
self.coated_substrate().astype(np.uint8),
cv2.RETR_EXTERNAL,
cv2.CHAIN_APPROX_SIMPLE,
)
vec = cnt - boundary_pts.transpose(1, 0, 2)
(I0, I1) = np.argmin(np.linalg.norm(vec, axis=-1), axis=0)
return np.squeeze(cnt[I0 : I1 + 1], axis=1)
[docs]
@attrcache("_average_thickness")
def average_thickness(self):
"""Average thickness of the coating layer.
Examples
--------
.. only:: doctest
>>> coat = getfixture('RectIfdRoughness_setup')
>>> coat.average_thickness()
50.25...
"""
idxs = self.interface_indices()
if len(idxs) == 0:
return np.nan
i0, i1 = idxs
subst_cnt = self.substrate_contour()[i0:i1]
A = np.count_nonzero(self.extract_layer())
(t,) = root(partial(_uniform_layer_area, subst=subst_cnt, x0=A), [0]).x
return float(t)
[docs]
def analyze(self):
"""Return analysis result.
Returns
-------
:class:`RectIfdRoughnessData`
"""
return self.DataType(self.average_thickness(), self.roughness()[0])
[docs]
def draw(self, pairs_dist=20.0):
"""Visualize the analysis result.
Draws the surface, the uniform layer, and the roughness pairs.
Parameters
----------
pairs_dist : float
Distance between the roughness pairs in the IFD parameter space.
Decreasing this value increases the density of pairs.
"""
image = cv2.cvtColor(self.image, cv2.COLOR_GRAY2RGB)
if not self.valid():
return image
image[self.extract_layer()] = 255
cv2.polylines(
image,
[self.surface().reshape(-1, 1, 2).astype(np.int32)],
isClosed=False,
color=(0, 0, 255),
thickness=1,
)
cv2.polylines(
image,
[self.uniform_layer().reshape(-1, 1, 2).astype(np.int32)],
isClosed=False,
color=(255, 0, 0),
thickness=1,
)
_, path = self.roughness()
path_len = np.sum(np.linalg.norm(np.diff(path, axis=0), axis=-1), axis=0)
u = np.linspace(0, path_len, int(path_len // pairs_dist))
pairs = sample_polyline(path, u)
pairs_surf_pt = sample_polyline(self.surface(), pairs[:, 0])
pairs_ul_pt = sample_polyline(self.uniform_layer(), pairs[:, 1])
pairs_pts = np.stack([pairs_surf_pt, pairs_ul_pt])[np.newaxis, ...]
cv2.polylines(
image,
pairs_pts.astype(np.int32).transpose(2, 1, 0, 3),
isClosed=False,
color=(0, 255, 0),
thickness=1,
)
return image
def _uniform_layer_area(thickness, subst, x0):
cnt = np.concatenate([subst, np.flip(parallel_curve(subst, thickness[0]), axis=0)])
return cv2.contourArea(cnt.astype(np.float32)) - x0