OpenDrop in a script

opendrop_script_example.py

from math import pi as PI

import cv2
import matplotlib.pyplot as plt

from opendrop.features.pendant import extract_pendant_features
from opendrop.fit.younglaplace import young_laplace_fit
from opendrop.geometry import Rect2


# Image to analyse
image = cv2.imread("/home/eugene/Projects/opendrop/example_images/water_in_air001.png")

# SI units
drop_density = 1000.0  # kg/m3
continuous_density = 0.0  # kg/m3
needle_diameter = 0.7176e-3  # m
gravity = 9.81  # m/s^2

# ROIs (comment to select using GUI)
#  needle_region = Rect2(x0=343, y0=34, x1=667, y1=136)
#  drop_region = Rect2(x0=214, y0=183, x1=790, y1=704)

# GUI options (comment to disable)
show_needle = True  # Show cropped needle image
show_drop = True  # Show cropped drop image
show_features = True  # Show extracted needle and drop edges
show_fit = True  # Show fitted drop profile and residuals


###############
# Select ROIs #
###############

if 'needle_region' not in globals():
    x, y, w, h = cv2.selectROI("Select needle region", image, fromCenter=False)
    # We use our own rect structure...
    needle_region = Rect2(x=x, y=y, w=w, h=h)
    cv2.destroyAllWindows()
    print("Needle region:", needle_region)
if 'drop_region' not in globals():
    x, y, w, h = cv2.selectROI("Select drop region", image, fromCenter=False)
    drop_region = Rect2(x=x, y=y, w=w, h=h) 
    cv2.destroyAllWindows()
    print("Drop region:", drop_region)

if globals().get('show_needle'):
    cv2.imshow(
        "Needle preview",
        image[
            needle_region.y0:needle_region.y1+1,
            needle_region.x0:needle_region.x1+1,
        ],
    )
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if globals().get('show_drop'):
    cv2.imshow(
        "Drop preview",
        image[
            drop_region.y0:drop_region.y1+1,
            drop_region.x0:drop_region.x1+1,
        ],
    )
    cv2.waitKey(0)
    cv2.destroyAllWindows()


####################
# Extract features #
####################
features = extract_pendant_features(
    image,
    drop_region=drop_region,
    needle_region=needle_region,
    # True to include mask of drop and needle pixels in the return value.
    labels=globals().get('show_features', False),
)

if globals().get('show_features'):
    labelled_image = image.copy()
    labelled_image[features.labels==1] = [255, 0, 0]  # Drop is blue
    labelled_image[features.labels==2] = [0, 255, 0]  # Needle is green
    cv2.imshow(
        "Extracted features",
        labelled_image,
    )
    cv2.waitKey(0)
    cv2.destroyAllWindows()


#############################
# Fit the extracted profile #
#############################
fit = young_laplace_fit(
    features.drop_points,
    verbose=True,
)

# Calculate physical quantities.
px_size = needle_diameter/features.needle_diameter
delta_density = abs(drop_density - continuous_density)
radius = fit.radius * px_size
surface_area = fit.surface_area * px_size**2
volume = fit.volume * px_size**3
ift = delta_density * gravity * radius**2 / fit.bond
worthington = (delta_density * gravity * volume) / (PI * ift * needle_diameter)

print(
    "\nResults\n=======\n",
    "Drop radius [mm]: ", radius*1e3, "\n",
    "Surface area [mm^2]: ", surface_area*1e6, "\n",
    "Volume [mm^3]: ", volume*1e9, "\n",
    "IFT [mN/m]: ", ift*1e3, "\n",
    "Worthington: ", worthington,
)


if globals().get('show_fit'):
    fig = plt.figure()
    ax1 = fig.add_subplot(1, 2, 1)
    ax1.set_title("Fitted drop profile")
    ax1.imshow(image)

    # fit.closest is a 2D array with shape (2, N) where N is the number of extracted points of the drop
    # profile.
    ax1.scatter(*fit.closest, s=0.5)

    ax2 = fig.add_subplot(1, 2, 2)
    ax2.set_title("Residuals")
    ax2.scatter(
        # Arclength from drop apex
        fit.arclengths,
        # Residual of fit, +ve means fitted profile is inside drop I think, and -ve means outside.
        fit.residuals,
        s=0.5,
    )

    plt.show(block=True)