Source code for pwspy.utility.machineVision
# Copyright 2018-2020 Nick Anthony, Backman Biophotonics Lab, Northwestern University
#
# This file is part of PWSpy.
#
# PWSpy is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PWSpy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with PWSpy. If not, see <https://www.gnu.org/licenses/>.
"""
Useful functions for processing images. Currently its contents are focused
on image stabilization.
Functions
------------
.. autosummary::
:toctree: generated/
to8bit
SIFTRegisterTransform
ORBRegisterTransform
edgeDetectRegisterTranslation
"""
from __future__ import annotations
import logging
import typing
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from skimage import feature
from skimage import registration
if typing.TYPE_CHECKING:
import cv2
logger = logging.getLogger(__name__)
[docs]def to8bit(arr: np.ndarray) -> np.ndarray:
"""Converts boolean or float type numpy arrays to 8bit and scales the data to span from 0 to 255. Used for many
OpenCV functions.
Args:
arr: The input array
Returns:
The output array of dtype numpy.uint8
"""
if arr.dtype == bool:
arr = arr.astype(np.uint8) * 255
else:
arr = arr.astype(float) # This solves problems if the array is int16 type
Min = np.percentile(arr, 0.1)
arr -= Min
Max = np.percentile(arr, 99.9)
arr = arr / Max * 255
arr[arr < 0] = 0
arr[arr > 255] = 255
return arr.astype(np.uint8)
[docs]def SIFTRegisterTransform(reference: np.ndarray, other: typing.Iterable[np.ndarray], mask: np.ndarray = None, debugPlots: bool = False) -> typing.Tuple[typing.List[np.ndarray], animation.ArtistAnimation]:
"""Given a 2D reference image and a list of other images of the same scene but shifted a bit this function will use OpenCV to calculate the transform from
each of the other images to the reference. The transforms can be inverted using cv2.invertAffineTransform().
It will return a list of transforms. Each transform is a 2x3 array in the form returned
by opencv.estimateAffinePartial2d(). a boolean mask can be used to select which areas will be searched for features to be used
in calculating the transform. This seems to work much better for normalized images.
This code is basically a copy of this example, it can probably be improved upon:
https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_feature_homography/py_feature_homography.html
Args:
reference (np.ndarray): The 2d reference image.
other (Iterable[np.ndarray]): An iterable containing the images that you want to calculate the translations for.
mask (np.ndarray): A boolean array indicating which parts of the reference image should be analyzed. If `None` then the whole image will be used.
debugPlots (bool): Indicates if extra plots should be openend showing the process of the function.
Returns:
tuple: A tuple containing:
List[np.ndarray]: Returns a list of transforms. Each transform is a 2x3 array in the form returned by opencv.estimateAffinePartial2d(). Note that even
though they are returned as affine transforms they will only contain translation information, no scaling, shear, or rotation.
ArtistAnimation: A reference the animation used to diplay the results of the function.
"""
import cv2
refImg = to8bit(reference)
if mask:
mask = mask.astype(np.uint8)
#Set up flann matcher
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv2.FlannBasedMatcher(index_params, search_params)
# Initiate SIFT detector
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(refImg, mask=mask)
transforms = []
if debugPlots:
anFig, anAx = plt.subplots()
anims = []
for i, img in enumerate(other):
logger.debug(f"Calculating SIFT matches for image {i} of {len(other)}")
otherImg = to8bit(img)
# find the keypoints and descriptors with SIFT
kp2, des2 = sift.detectAndCompute(otherImg, mask=None)
good = _knnMatch(flann, des1, des2)
MIN_MATCH_COUNT = 5
if len(good) > MIN_MATCH_COUNT:
src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)
M, inlierMask = cv2.estimateAffinePartial2D(src_pts, dst_pts)
matchesMask = inlierMask.ravel().tolist()
else:
logger.warning(f"Image {i}: Not enough matches are found - {len(good)}/{MIN_MATCH_COUNT}")
M = None
matchesMask = None
transforms.append(M)
if debugPlots and (M is not None):
anims.append([anAx.imshow(cv2.warpAffine(otherImg, cv2.invertAffineTransform(M), otherImg.shape), 'gray')])
h, w = refImg.shape
pts = np.float32([[0, 0], [0, h - 1], [w - 1, h - 1], [w - 1, 0]]).reshape(-1, 1, 2)
dst = cv2.transform(pts, M)
draw_params = dict(matchColor=(0, 255, 0), # draw matches in green color
singlePointColor=None,
matchesMask=matchesMask, # draw only inliers
flags=2)
img2 = cv2.polylines(otherImg, [np.int32(dst)], True, 255, 3, cv2.LINE_AA)
img3 = cv2.drawMatches(refImg, kp1, img2, kp2, good, None, **draw_params)
fig, ax = plt.subplots()
ax.imshow(img3, 'gray')
if debugPlots:
anFig.suptitle("SIFT: If transforms worked, image should not appear to move.")
an = animation.ArtistAnimation(anFig, anims)
else:
an = None
return transforms, an
[docs]def ORBRegisterTransform(reference: np.ndarray, other: typing.Iterable[np.ndarray], mask: np.ndarray = None, debugPlots: bool = False) -> typing.Tuple[typing.List[np.ndarray], animation.ArtistAnimation]:
"""Given a 2D reference image and a list of other images of the same scene but shifted a bit this function will use OpenCV to calculate the transform from
each of the other images to the reference. The transforms can be inverted using cv2.invertAffineTransform().
It will return a list of transforms. Each transform is a 2x3 array in the form returned
by opencv.estimateAffinePartial2d(). a boolean mask can be used to select which areas will be searched for features to be used
in calculating the transform.
Args:
reference (np.ndarray): The 2d reference image.
other (Iterable[np.ndarray]): An iterable containing the images that you want to calculate the translations for.
mask (np.ndarray): A boolean array indicating which parts of the reference image should be analyzed. If `None` then the whole image will be used.
debugPlots (bool): Indicates if extra plots should be openend showing the process of the function.
Returns:
tuple: A tuple containing:
List[np.ndarray]: Returns a list of transforms. Each transform is a 2x3 array in the form returned by opencv.estimateAffinePartial2d(). Note that even
though they are returned as affine transforms they will only contain translation information, no scaling, shear, or rotation.
ArtistAnimation: A reference the animation used to diplay the results of the function.
"""
import cv2
refImg = to8bit(reference)
if mask:
mask = mask.astype(np.uint8)
#Create FLANN matcher
FLANN_INDEX_LSH = 6
index_params = dict(algorithm=FLANN_INDEX_LSH,
table_number=6, # 12
key_size=12, # 20
multi_probe_level=1) # 2
search_params = dict(checks=100)
flann = cv2.FlannBasedMatcher(index_params, search_params)
# Initiate ORB detector
orb = cv2.ORB_create()
kp1, des1 = orb.detectAndCompute(refImg, mask=mask)
transforms = []
if debugPlots:
anFig, anAx = plt.subplots()
anims = []
for i, img in enumerate(other):
otherImg = to8bit(img)
# find the keypoints and descriptors with ORB
kp2, des2 = orb.detectAndCompute(otherImg, mask=None)
good = _knnMatch(flann, des1, des2)
MIN_MATCH_COUNT = 5
if len(good) > MIN_MATCH_COUNT:
src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)
M, inlierMask = cv2.estimateAffinePartial2D(src_pts, dst_pts)
matchesMask = inlierMask.ravel().tolist()
else:
logging.getLogger(__name__).warning(f"Image {i}: Not enough matches are found - {len(good)}/{MIN_MATCH_COUNT}")
M = None
matchesMask = None
transforms.append(M)
if debugPlots and (M is not None):
anims.append([anAx.imshow(cv2.warpAffine(otherImg, cv2.invertAffineTransform(M), otherImg.shape), 'gray')])
h, w = refImg.shape
pts = np.float32([[0, 0], [0, h - 1], [w - 1, h - 1], [w - 1, 0]]).reshape(-1, 1, 2)
dst = cv2.transform(pts, M)
draw_params = dict(matchColor=(0, 255, 0), # draw matches in green color
singlePointColor=None,
matchesMask=matchesMask, # draw only inliers
flags=2)
img2 = cv2.polylines(otherImg, [np.int32(dst)], True, 255, 3, cv2.LINE_AA)
img3 = cv2.drawMatches(refImg, kp1, img2, kp2, good, None, **draw_params)
fig, ax = plt.subplots()
ax.imshow(img3, 'gray')
if debugPlots:
anFig.suptitle("ORB: If transforms worked, image should not appear to move.")
an = animation.ArtistAnimation(anFig, anims)
else:
an = None
return transforms, an
def _knnMatch(flannMatcher: cv2.FlannBasedMatcher, des1: np.ndarray, des2: np.ndarray) -> typing.List[cv2.DMatch]:
"""
Return a list of sufficiently good matches between keypoint descriptors.
Args:
flannMatcher: The opencv FLANN matcher object
des1: An array of `query descriptors` to find matches for.
des2: An array of `train descriptors` to check as matches for des1
Returns:
A list of the selecte opencv DMatch objects
"""
matches = flannMatcher.knnMatch(des1, des2, k=2) # Return the two best matches. Note, might only return 1 match
# store all the good matches as per Lowe's ratio test.
good = []
for iMatches in matches:
if len(iMatches) == 1:
good.append(iMatches[0]) # If only one match was found keep it
elif len(iMatches) == 2:
if iMatches[0].distance < 0.7 * iMatches[1].distance: # This is known as Lowe's Ratio Test. If two matches were returned only keep the best match if it is sufficiently better than the second best match
good.append(iMatches[0])
else:
raise Exception("Programming Error!")
return good
[docs]def edgeDetectRegisterTranslation(reference: np.ndarray, other: typing.Iterable[np.ndarray], mask: np.ndarray = None, debugPlots: bool = False, sigma: float = 3) -> typing.Tuple[typing.Iterable[np.ndarray], typing.List]:
"""This function is used to find the relative translation between a reference image and a list of other similar images. Unlike `SIFRegisterTransforms` this function
will not work for images that are rotated relative to the reference. However, it does provide more robust performance for images that do not look identical.
Args:
reference (np.ndarray): The 2d reference image.
other (Iterable[np.ndarray]): An iterable containing the images that you want to calculate the translations for.
mask (np.ndarray): A boolean array indicating which parts of the reference image should be analyzed. If `None` then the whole image will be used.
debugPlots (bool): Indicates if extra plots should be openend showing the process of the function.
sigma (float): this parameter is passed to skimage.feature.canny to detect edges.
Returns:
tuple: A tuple containing:
list[np.ndarray]: Returns a list of transforms. Each transform is a 2x3 array in the form returned by opencv.estimateAffinePartial2d(). Note that even
though they are returned as affine transforms they will only contain translation information, no scaling, shear, or rotation.
list: A list of references to the plotting widgets used to display the results of the function.
"""
import cv2
from mpl_qt_viz.visualizers import MultiPlot
refEd = feature.canny(reference, sigma=sigma)
if mask is not None: refEd[~mask] = False # Clear any detected edges outside of the mask
imEd = [feature.canny(im, sigma=sigma) for im in other]
affineTransforms = []
if debugPlots:
anEdFig, anEdAx = plt.subplots()
anEdFig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
anEdAx.get_xaxis().set_visible(False)
anEdAx.get_yaxis().set_visible(False)
anFig, anAx = plt.subplots()
anFig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
anAx.get_xaxis().set_visible(False)
anAx.get_yaxis().set_visible(False)
anims = [[anAx.imshow(to8bit(reference), 'gray'), anAx.text(100, 100, "Reference", color='r')]]
animsEd = [[anEdAx.imshow(to8bit(refEd), 'gray'), anEdAx.text(100, 100, "Reference", color='w')]]
for i, (im, edgeIm) in enumerate(zip(other, imEd)):
shifts, error, phasediff = feature.register_translation(edgeIm, refEd)
logging.getLogger(__name__).info(f"Translation: {shifts}, RMS Error: {error}, Phase Difference:{phasediff}")
shifts = np.array([[1, 0, shifts[1]],
[0, 1, shifts[0]]], dtype=float) # Convert the shift to an affine transform
affineTransforms.append(shifts)
if debugPlots:
plt.figure()
animsEd.append([anEdAx.imshow(cv2.warpAffine(to8bit(edgeIm), cv2.invertAffineTransform(shifts), edgeIm.shape), 'gray'), anEdAx.text(100, 100, str(i), color='w')])
anims.append([anAx.imshow(cv2.warpAffine(to8bit(im), cv2.invertAffineTransform(shifts), im.shape), 'gray'), anAx.text(100, 100, str(i), color='r')])
if debugPlots:
an = [MultiPlot(anims, "If transforms worked, cells should not appear to move."), MultiPlot(animsEd, "If transforms worked, cells should not appear to move.")]
[i.show() for i in an]
else:
an = []
return affineTransforms, an
def crossCorrelateRegisterTranslation(reference: np.ndarray, other: typing.Iterable[np.ndarray], debugPlots: bool = False) -> typing.Tuple[typing.Iterable[np.ndarray], 'mpl_qt_viz.visualizers.MultiPlot']:
"""This function is used to find the relative translation between a reference image and a list of other similar images. Unlike `SIFRegisterTransforms` this function
will not work for images that are rotated relative to the reference.
Args:
reference (np.ndarray): The 2d reference image.
other (Iterable[np.ndarray]): An iterable containing the images that you want to calculate the translations for.
debugPlots (bool): Indicates if extra plots should be openend showing the process of the function.
Returns:
tuple: A tuple containing:
list[np.ndarray]: Returns a list of transforms. Each transform is a 2x3 array in the form returned by opencv.estimateAffinePartial2d(). Note that even
though they are returned as affine transforms they will only contain translation information, no scaling, shear, or rotation.
MultiPlot: A reference to the plotting widgets used to display the results of the function. If `debugPlots` is False this will be `None`
"""
import cv2
if debugPlots:
anFig, anAx = plt.subplots()
anFig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
anAx.get_xaxis().set_visible(False)
anAx.get_yaxis().set_visible(False)
anims = [[anAx.imshow(to8bit(reference), 'gray'), anAx.text(100, 100, "Reference", color='r')]]
affineTransforms = []
for i, im in enumerate(other):
shifts, error, phasediff = registration.phase_cross_correlation(im, reference, return_error=True)
logging.getLogger(__name__).debug(f"Translation: {shifts}, RMS Error: {error}, Phase Difference:{phasediff}")
shifts = np.array([[1, 0, shifts[1]],
[0, 1, shifts[0]]], dtype=float) # Convert the shift to an affine transform
affineTransforms.append(shifts)
if debugPlots:
anims.append([
anAx.imshow(cv2.warpAffine(to8bit(im), cv2.invertAffineTransform(shifts), im.shape), 'gray'),
anAx.text(100, 100, str(i), color='r')])
if debugPlots:
from mpl_qt_viz.visualizers import MultiPlot
an = MultiPlot(anims, "If transforms worked, images should not appear to move.")
an.show()
else:
an = None
return affineTransforms, an