rahix/openflexure-microscope-server
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

First working camera stage mapping calibration

Browse source
This commit is contained in:
Richard Bowman 2020-03-12 16:03:11 +00:00
parent c240453762
commit e452d623ca
8 changed files with 1078 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

1
openflexure_microscope/api/default_extensions/__init__.py
View file

@ -5,6 +5,7 @@ from .autofocus import autofocus_extension_v2
from .scan import scan_extension_v2
from .zip_builder import zip_extension_v2
from .autostorage import autostorage_extension_v2
from .camera_stage_mapping import csm_extension
# "Gracefully" handle cases where picamera cannot be imported (eg test server)
try:

1
openflexure_microscope/api/default_extensions/camera_stage_mapping/__init__.py Normal file
View file

@ -0,0 +1 @@
from .extension import csm_extension

95
openflexure_microscope/api/default_extensions/camera_stage_mapping/array_with_attrs.py Normal file
View file

@ -0,0 +1,95 @@
# -*- coding: utf-8 -*-
"""
Created on Tue May 26 08:08:14 2015
@author: rwb27
"""
import numpy as np
class AttributeDict(dict):
"""This class extends a dictionary to have a "create" method for
compatibility with h5py attrs objects."""
def create(self, name, data):
self[name] = data
def modify(self, name, data):
self[name] = data
def copy_arrays(self):
"""Replace any numpy.ndarray in the dict with a copy, to break any unintentional links."""
for k in list(self.keys()):
if isinstance(self[k], np.ndarray):
self[k] = np.copy(self[k])
def ensure_attribute_dict(obj, copy=False):
"""Given a mapping that may or not be an AttributeDict, return an
AttributeDict object that either is, or copies the data of, the input."""
if isinstance(obj, AttributeDict) and not copy:
return obj
else:
out = AttributeDict(obj)
if copy:
out.copy_arrays()
return out
def ensure_attrs(obj):
"""Return an ArrayWithAttrs version of an array-like object, may be the
original object if it already has attrs."""
if hasattr(obj, 'attrs'):
return obj #if it has attrs, do nothing
else:
return ArrayWithAttrs(obj) #otherwise, wrap it
class ArrayWithAttrs(np.ndarray):
"""A numpy ndarray, with an AttributeDict accessible as array.attrs.
This class is intended as a temporary version of an h5py dataset to allow
the easy passing of metadata/attributes around nplab functions. It owes
a lot to the ``InfoArray`` example in `numpy` documentation on subclassing
`numpy.ndarray`.
"""
def __new__(cls, input_array, attrs={}):
"""Make a new ndarray, based on an existing one, with an attrs dict.
This function adds an attributes dictionary to a numpy array, to make
it work like an h5py dataset. It doesn't copy data if it can be
avoided."""
# the input array should be a numpy array, then we cast it to this type
obj = np.asarray(input_array).view(cls)
# next, add the dict
# ensure_attribute_dict always returns an AttributeDict
obj.attrs = ensure_attribute_dict(attrs)
# return the new object
return obj
def __array_finalize__(self, obj):
# this is called by numpy when the object is created (__new__ may or
# may not get called)
if obj is None: return # if obj is None, __new__ was called - do nothing
# if we didn't create the object with __new__, we must add the attrs
# dictionary. We copy this from the source object if possible (while
# ensuring it's the right type) or create a new, empty one if not.
# NB we don't use ensure_attribute_dict because we want to make sure the
# dict object is *copied* not merely referenced.
self.attrs = ensure_attribute_dict(getattr(obj, 'attrs', {}), copy=True)
def attribute_bundler(attrs):
"""Return a function that bundles the supplied attributes with an array."""
def bundle_attrs(array):
return ArrayWithAttrs(array, attrs=attrs)
class DummyHDF5Group(dict):
def __init__(self,dictionary, attrs ={}, name="DummyHDF5Group"):
super(DummyHDF5Group, self).__init__()
self.attrs = attrs
for key in dictionary:
self[key] = dictionary[key]
self.name = name
self.basename = name
file = None
parent = None

279
openflexure_microscope/api/default_extensions/camera_stage_mapping/camera_stage_calibration_1d.py Normal file
View file

@ -0,0 +1,279 @@
"""
1D calibration of the relationship between a stage and a camera
The `Tracker` class in this file is used to simplify code for tasks that involve moving the
stage, and tracking the corresponding motion with the camera.
(c) Richard Bowman 2019, released under GNU GPL v3
"""
import numpy as np
import time
from numpy.linalg import norm
from .camera_stage_tracker import Tracker, move_until_motion_detected
import logging
def displacements(positions):
"""Calculate the absolute distance of each point from the first point."""
return norm(positions - positions[0,:][np.newaxis,:], axis=1)
def direction_from_points(points):
"""Given an Nx2 array of points, figure out the principal component.
The return value is a normalised vector that points along the
direction with the most motion.
"""
points = points.astype(np.float)
points -= np.mean(points, axis=0)[np.newaxis, :]
eigenvalues, eigenvectors = np.linalg.eig(np.cov(points.T))
return eigenvectors[:,np.argmax(eigenvalues)]
def apply_backlash(x, backlash=0, start_unwound=True):
"""Apply a basic model of backlash to a set of coordinates.
The output (y) will lag behind the input by up to `backlash`
`start_unwound` (default: True) assumes we change direction
at the start of the time series, so you will get no motion
until `x[i]` has moved by at least `2*backlash`.
"""
y = np.zeros_like(x)
if start_unwound:
initial_direction = np.sign(x[1] - x[0])
y[0] = x[0] + initial_direction * backlash
else:
y[0] = x[0]
for i in range(1,len(x)):
d = x[i] - y[i-1]
if np.abs(d) >= backlash:
y[i] = x[i] - np.sign(d) * backlash
else:
y[i] = y[i-1]
return y
def fit_backlash(moves):
"""Given a set of linear moves forwards and back, estimate backlash.
The result is an estimate of the amount of backlash, and the ratio
of steps to pixels. The moves should be in the same
format as `Tracker.history`.
We use a very basic fitting method: we do a brute-force search for
the backlash value, and for each value of backlash we fit a line to
the relationship between stage position (after modelling backlash)
and image position. We then pick the value of backlash that gets
the lowest residuals. Currently the backlash values tried will
start at 0 and increase by 1 or by a factor of 1.33 each time.
The return value is a dictionary with the following keys:
backlash: float
the estimated backlash, in motor steps
pixels_per_step: float
the gradient of pixels to steps
fractional_error: float
an estimate of the goodness of fit
stage_direction: numpy.ndarray
unit vector in the direction of stage motion
image_direction: numpy.ndarray
unit vector in the direction of the motion measured
on the camera
pixels_per_step_vector: numpy.ndarray
The displacement in 2D on the camera resulting from
one step in `stage_direction`. This is equal to the
product of `pixels_per_step` and `image_direction`.
"""
all_stage_points, all_image_points = moves
# Figure out the direction of motion, and reduce everything to 1D
image_direction = direction_from_points(all_image_points)
stage_direction = direction_from_points(all_stage_points)
xfit = np.sum(all_stage_points * stage_direction[np.newaxis, :], axis=1)
yfit = np.sum(all_image_points * image_direction[np.newaxis, :], axis=1)
# We should probably use a fancy optimiser to fit the backlash, but
# brute-forcing it is reliable and doesn't take long.
def fit_motion(xfit, yfit, backlash=0):
"""Using the model of backlash, fit the observed camera motion"""
xfit_blsh = apply_backlash(xfit, backlash)
xfit_blsh -= np.mean(xfit_blsh)
m, c = np.polyfit(xfit_blsh, yfit, 1)
residuals = yfit - (xfit_blsh * m + c)
return m, c, np.std(residuals, ddof=3)
max_backlash = (np.max(xfit) - np.min(xfit))/3
backlash_values = []
residual_values = []
backlash = 0
while backlash < max_backlash:
m, c, residual = fit_motion(xfit, yfit, backlash)
residual_values.append(residual)
backlash_values.append(backlash)
backlash += max(1, backlash/3)
backlash = backlash_values[np.argmin(residual_values)]
m, c, residual = fit_motion(xfit, yfit, backlash)
fractional_error = residual/norm(np.diff(yfit))
if fractional_error > 0.1:
raise ValueError("The fit didn't look successful")
return {
"backlash": backlash,
"pixels_per_step": m,
"fractional_error": fractional_error,
"stage_direction": stage_direction,
"image_direction": image_direction,
"pixels_per_step_vector": m * image_direction,
}
def calibrate_backlash_1d(tracker, move, direction=np.array([1,0,0])):
"""Figure out reasonable step sizes for calibration, and estimate the backlash."""
try: # Ensure that the tracker has a template set
_ = tracker.template
except:
tracker.acquire_template()
assert tracker.stage_positions.shape[0] == 1
original_stage_pos = tracker.stage_positions[-1,:]
direction = direction / np.sum(direction**2)**0.5 # ensure "direction" is normalised
logging.info("Moving the stage until we see motion...")
# Move the stage until we can see a significant amount of motion
i, m = move_until_motion_detected(
tracker, move, direction, threshold=tracker.max_safe_displacement * 0.2)
logging.info("Moving the stage to the edge of the field of view...")
i, m = move_until_motion_detected(
tracker, move, direction,
threshold=tracker.max_safe_displacement * 0.7,
multipliers=m/2.0 * np.arange(20),
detect_cumulative_motion=True)
exponential_moves = tracker.history
# Include this final step, and make a rough estimate of the scaling from stage to image
stage_pos, image_pos = tracker.history
stage_step = stage_pos[-1, :] - stage_pos[-1 - i, :]
image_step = image_pos[-1, :] - image_pos[-1 - i, :]
steps_per_pixel = norm(stage_step)/norm(image_step)
# Calculate a step that moves roughly 0.2 times the max. displacement (i.e. 0.1 times the FoV)
sensible_step = direction * tracker.max_safe_displacement * 0.2 * steps_per_pixel
tracker.reset_history()
logging.info("Moving the stage backwards to measure backlash (1/2)")
# Now move backwards, in 10 steps that should roughly cross the field of view.
# If the stage has no backlash, this will move too far, hence the break statement to
# prevent it moving outside of the field of view.
starting_stage_pos, starting_camera_pos = tracker.append_point()
for i in range(15):
move(starting_stage_pos - sensible_step * (i + 1))
#print(".", end="")
stage_pos, image_pos = tracker.append_point()
if (i > 3 and tracker.moving_away_from_centre
and norm(image_pos) > 0.65 * tracker.max_safe_displacement):
break # Stop once we have moved far enough
logging.info("Moving the stage forwards to measure backlash (2/2)")
# Move forwards again, in 10 steps
starting_stage_pos, starting_camera_pos = tracker.append_point()
for i in range(15):
move(starting_stage_pos + sensible_step * (i + 1))
#print(".", end="")
stage_pos, image_pos = tracker.append_point()
if (i > 3 and tracker.moving_away_from_centre
and norm(image_pos) > 0.65 * tracker.max_safe_displacement):
break # Stop once we have moved far enough
linear_moves = tracker.history
try:
res = fit_backlash(linear_moves)
backlash_correction = sensible_step / norm(sensible_step) * res["backlash"] * 1.5
# Finally, move back to the starting position, doing backlash-corrected moves.
logging.info("Moving back to the start, correcting for backlash...")
tracker.reset_history()
stage_pos, camera_pos = tracker.append_point()
while np.dot(stage_pos - sensible_step - original_stage_pos, sensible_step) > 0:
move(stage_pos - sensible_step - backlash_correction)
move(stage_pos - sensible_step)
stage_pos, camera_pos = tracker.append_point()
backlash_corrected_moves = tracker.history
move(original_stage_pos - backlash_correction)
except ValueError:
return {"exponential_moves": exponential_moves, "linear_moves": linear_moves,}
finally:
# Reset position
move(original_stage_pos)
logging.info(f"Estimated backlash {res['backlash']:.0f} steps")
logging.info(f"Stage-to-image ratio {np.abs(res['pixels_per_step']):.3f} pixels/step")
logging.info(f"Residuals were about {res['fractional_error']:.2f} times the step size")
res.update({
"exponential_moves": exponential_moves,
"linear_moves": linear_moves,
"backlash_corrected_moves": backlash_corrected_moves
})
return res
def plot_1d_backlash_calibration(results):
"""Plot the results of a calibration run"""
from matplotlib import pyplot as plt
f, ax = plt.subplots(1,2)
for k in ["exponential", "linear", "backlash_corrected"]:
moves = results[k+"_moves"]
if moves is not None:
ax[0].plot(moves[1][:,0], moves[1][:,1], 'o-')
ax[0].set_aspect(1, adjustable="datalim")
image_direction = results["image_direction"]
stage_direction = results["stage_direction"]
def convert_moves(moves):
stage_pos, image_pos = moves
stage_1d = np.sum(stage_pos * stage_direction[np.newaxis, :], axis=1)
image_1d = np.sum(image_pos * image_direction[np.newaxis, :], axis=1)
return stage_1d, image_1d
ax[1].plot(*convert_moves(results["exponential_moves"]), 'o-')
stage_pos, image_pos = convert_moves(results["linear_moves"])
model = apply_backlash(stage_pos, results["backlash"])
model *= results["pixels_per_step"]
model += np.mean(image_pos) - np.mean(model)
ax[1].plot(stage_pos, model, '-')
ax[1].plot(stage_pos, image_pos, 'o')
if results["backlash_corrected_moves"] is not None:
ax[1].plot(*convert_moves(results["backlash_corrected_moves"]), '+')
return f, ax
def image_to_stage_displacement_from_1d(calibrations):
"""Combine X and Y calibrations
This uses the output from `calibrate_backlash_1d`, run at least
twice with orthogonal (or at least different) `direction` parameters.
The resulting 2x2 transformation matrix should map from image
to stage coordinates. Currently, the backlash estimate given
by this function is only really trustworthy if you've supplied
two orthogonal calibrations - that will usually be the case.
"""
stage_vectors = []
image_vectors = []
backlash = np.zeros(3)
for cal in calibrations:
stage_vectors.append(cal["stage_direction"][:2])
image_vectors.append(cal["pixels_per_step_vector"])
# our backlash estimate will be the maximum backlash
# measured in each direction
c_blash = np.abs(cal["backlash"] * cal["stage_direction"])
backlash[backlash < c_blash] = c_blash[backlash < c_blash]
A, res, rank, s = np.linalg.lstsq(image_vectors, stage_vectors) # we solve image*A = stage
return {
"image_to_stage_displacement": A,
"backlash_vector": backlash,
"backlash": np.max(backlash),
}

81
openflexure_microscope/api/default_extensions/camera_stage_mapping/camera_stage_calibration_2d.py Normal file
View file

@ -0,0 +1,81 @@
"""
Camera-stage calibration, 2D
Uses 2D motion to try to calibrate the relationship between a camera and a stage.
(c) Richard Bowman 2019, released under GNU GPL v3
"""
import numpy as np
import time
from numpy.linalg import norm
from matplotlib import pyplot as plt
from camera_stage_tracker import Tracker, move_until_motion_detected
from functools import partial
def backlash_corrected_move(get_position, move, backlash_amount, pos):
"""Make two moves, arriving at `pos` from a consistent direction"""
displacement = pos - get_position()
backlash_vector = (displacement < 0).astype(np.int)*backlash_amount
if np.any(backlash_vector > 0):
move(pos - backlash_vector)
move(pos)
def bake_backlash_corrected_move(get_position, move, backlash_amount):
"""Return a function that performs backlash-corrected moves"""
return partial(backlash_corrected_move, get_position, move, backlash_amount)
def calibrate_xy_grid(tracker, move, step = 100, n_steps=4, backlash_compensation=0):
"""Make a series of moves in X and Y to determine the XY components of the pixel-to-sample matrix.
Arguments:
tracker: Tracker
An initialised Tracker object, centred on the starting point. This provides position readout from the stage and the camera.
move: function
A function that accepts a 1D array and performs an absolute move to
that position. If backlash correction is needed, include it here.
step : float, optional (default 100)
The amount to move the stage by. This should move the sample by approximately 1/10th of the field of view.
"""
try: # Ensure that the tracker has a template set
_ = tracker.template
except:
tracker.acquire_template()
tracker.reset_history() # make sure we get rid of the initial (0,0) point
starting_position = tracker.get_position()
# Move the stage in a square, recording the displacement from both the stage and the camera
try:
for x in (np.arange(n_steps) - n_steps/2.0)*step:
for y in (np.arange(n_steps) - n_steps/2.0)*step:
move(starting_position + np.array([x, y, 0]))
tracker.append_point()
finally:
move(starting_position)
# We then use least-squares to fit the XY part of the matrix relating
# pixels to distance
# stage_positions should be the stage positions, with a zero mean.
# image_positions should be the same, but calculated from the images
stage_positions, image_positions = tracker.history
stage_positions = stage_positions.astype(np.float)
stage_positions -= np.mean(stage_positions, axis=0)
stage_positions = stage_positions[:,:2] # ensure it's 2d
image_positions -= np.mean(image_positions, axis=0)
#image_positions *= -1 # To get the matrix right, we want the position of each
# image relative to the template, rather than the other way around
A, res, rank, s = np.linalg.lstsq(image_positions, stage_positions) # we solve pixel_shifts*A = location_shifts
transformed_image_positions = np.dot(image_positions, A)
residuals = transformed_image_positions - stage_positions
fractional_error = norm(residuals) / stage_positions.shape[0] step
print(f"Ratio of residuals to displacement is {fractional_error})")
if fractional_error > 0.05: # Check it was a reasonably good fit
print("Warning: the error fitting measured displacements was %.1f%%" % (fractional_error*100))
print(f"Calibrated the pixel-location matrix.\nResiduals were {fractional_error*100:.1f}% of the shift.")
return {
"image_to_stage_displacement": A,
"moves": (stage_positions, image_positions),
"fractional_error": fractional_error
}

290
openflexure_microscope/api/default_extensions/camera_stage_mapping/camera_stage_tracker.py Normal file
View file

@ -0,0 +1,290 @@
"""
Camera-stage tracker
The `Tracker` class in this file is used to simplify code for tasks that involve moving the
stage, and tracking the corresponding motion with the camera.
(c) Richard Bowman 2019, released under GNU GPL v3
"""
import numpy as np
import time
from numpy.linalg import norm
from matplotlib import pyplot as plt
import cv2
from scipy import ndimage
def central_half(image):
"""Return the central 50% (in X and Y) of an image"""
w, h = image.shape[:2]
return image[int(w/4):int(3*w/4),int(h/4):int(3*h/4), ...]
def datum_pixel(image):
"""Get the datum pixel of an image - if no property is present, assume the central pixel."""
try:
return np.array(image.datum_pixel)
except:
return (np.array(image.shape[:2]) - 1) / 2.
def locate_feature_in_image(image, feature, margin=0, restrict=False):
"""Find the given feature (small image) and return the position of its datum (or centre) in the image's pixels.
image : numpy.array
The image in which to look.
feature : numpy.array
The feature to look for. Ideally should be an `ImageWithLocation`.
margin : int (optional)
Make sure the feature image is at least this much smaller than the big image. NB this will take account of the
image datum points - if the datum points are superimposed, there must be at least margin pixels on each side of
the feature image.
restrict : bool (optional, default False)
If set to true, restrict the search area to a square of (margin * 2 + 1) pixels centred on the pixel that most
closely overlaps the datum points of the two images.
The `image` must be larger than `feature` by a margin big enough to produce a meaningful search area. We use the
OpenCV `matchTemplate` method to find the feature. The returned position is the position, relative to the corner of
the first image, of the "datum pixel" of the feature image. If no datum pixel is specified, we assume it's the
centre of the image. The output of this function can be passed into the pixel_to_location() method of the larger
image to yield the position in the sample of the feature you're looking for.
"""
# The line below is superfluous if we keep the datum-aware code below it.
assert image.shape[0] > feature.shape[0] and image.shape[1] > feature.shape[1], "Image must be larger than feature!"
# Check that there's enough space around the feature image
lower_margin = datum_pixel(image) - datum_pixel(feature)
upper_margin = (image.shape[:2] - datum_pixel(image)) - (feature.shape[:2] - datum_pixel(feature))
assert np.all(np.array([lower_margin, upper_margin]) >= margin), "The feature image is too large."
#TODO: sensible auto-crop of the template if it's too large?
image_shift = np.array((0,0))
if restrict:
# if requested, crop the larger image so that our search area is (2*margin + 1) square.
image_shift = np.array(lower_margin - margin,dtype = int)
image = image[image_shift[0]:image_shift[0] + feature.shape[0] + 2 * margin + 1,
image_shift[1]:image_shift[1] + feature.shape[1] + 2 * margin + 1, ...]
corr = cv2.matchTemplate(image, feature,
cv2.TM_SQDIFF_NORMED) # correlate them: NB the match position is the MINIMUM
corr = -corr # invert the image so we can find a peak
corr += (corr.max() - corr.min()) * 0.1 - corr.max() # background-subtract 90% of maximum
corr = cv2.threshold(corr, 0, 0, cv2.THRESH_TOZERO)[
1] # zero out any negative pixels - but there should always be > 0 nonzero pixels
assert np.sum(corr) > 0, "Error: the correlation image doesn't have any nonzero pixels."
peak = ndimage.measurements.center_of_mass(corr) # take the centroid (NB this is of grayscale values, not binary)
pos = np.array(peak) + image_shift + datum_pixel(feature) # return the position of the feature's datum point.
return pos
class Tracker():
def __init__(self, grab_image, get_position, settle=None):
"""A class to manage moving the stage and following motion in the image
Constructor Arguments:
grab_image: a function that returns the image as a numpy array
get_position: a function that returns position as a numpy array
settle: a function that waits and/or discards images
We accept functions because that seems like the easiest way to be
compatible with many different cameras/stages. Subclass and override
``__init__`` if you want to use a particular object instead.
NB the ``image_position`` that this class returns may be the negative of
what you might expect. This is because normally we are looking for
where a certain object (usually matched to a template image) is within
an image. Instead, the ``Tracker`` is following motion of the image
relative to a template. Our model is that we have a static image on
the slide, so we're tracking the slide's motion.
"""
self._grab_image = grab_image
self._get_position = get_position
self._settle = settle
self._template = None
self.margin = np.array([0, 0])
self._template_position = np.array([0.0, 0.0])
self.image_shape = None
def get_position(self):
"""Get the position of the stage"""
return np.array(self._get_position())
def settle(self):
"""Wait a short time and discard an image so the stage is no longer wobbling."""
if self._settle is not None:
self._settle()
else:
time.sleep(0.3)
self._grab_image()
@property
def template(self):
"""The template image (should be a numpy array)"""
if self._template is None:
raise ValueError("Attempt to use the tracker before setting the template")
else:
return self._template
@template.setter
def template(self, new_value):
self._template = new_value
def acquire_template(self, settle=True, reset_history=True, relative_positions=True):
"""Take a new image, and use it as the template. NB this records the initial point.
We will wait for the stage to settle, then acquire a new image to use as the template.
Immediately afterwards, we record the first point, so we will acquire a second image
and also read the stage's position.
The template image will be the central 50% of the starting image, which means the
maximum displacement will be 0.25 fields-of-view in all directions.
Arguments:
settle: bool, default True
Whether to wait for the stage to settle before taking the template image
reset_history: bool, default True
Whether to erase all the previously-stored positions
relative_positions: bool, default True
If true, we will define the first point (as read from the camera) to be [0,0]
and make all future measurements relative to this one. NB this won't affect
the stage positions, which are always absolute.
"""
if settle:
self.settle()
image = self._grab_image()
self.template = central_half(image)
self.image_shape = image.shape
self.margin = np.array(image.shape)[:2] - np.array(self.template.shape)[:2]
if reset_history:
self.reset_history()
self._template_position = np.array([0., 0.])
if relative_positions:
self._template_position = self.track_image(image) # Position should be zero initially
self.append_point(settle=False)
@property
def max_displacement(self):
"""The highest position values that can be tracked"""
return self.margin // 2 # TODO: be cleverer about non-trivial values of template_position
@property
def min_displacement(self):
"""The lowest position values that can be tracked"""
return -self.max_displacement # TODO: be cleverer about non-trivial template_position values
@property
def max_safe_displacement(self):
"""The biggest displacement we can safely attempt to track without knowing direction."""
return np.min(np.concatenate([self.max_displacement, -self.min_displacement]))
def track_image(self, image):
"""Find the position of the image relative to the template
NB this class is intended to track motion of the sample - most of
the time, we're interested in the motion of a (small) object that
is represented by the template, relative to the (larger) image. In
our case, we're doing the opposite - tracking motion of the image,
relative to a picture of part of the sample. That's why there is
a minus sign in front of `locate_feature_in_image` in the source
code.
"""
return - locate_feature_in_image(image, self.template) - self._template_position
def append_point(self, settle=True, image=None):
"""Find the current position using both stage and image, and append it"""
if settle:
self.settle()
if image is None:
image = self._grab_image()
image_pos = self.track_image(image)
stage_pos = self.get_position()
self._image_positions.append(image_pos)
self._stage_positions.append(stage_pos)
return stage_pos, image_pos
@property
def stage_positions(self):
"""An array of positions we have moved the stage to"""
return np.array(self._stage_positions)
@property
def image_positions(self):
"""An array of positions we have moved the stage to"""
return np.array(self._image_positions)
@property
def history(self):
"""Return arrays of stage, image positions"""
return self.stage_positions, self.image_positions
def reset_history(self, leave_first_point=False):
"""Reset the positions and displacements recorded"""
if leave_first_point:
self._stage_positions = [self._stage_positions[0]]
self._image_positions = [self._image_positions[0]]
else:
self._stage_positions = []
self._image_positions = []
@property
def moving_away_from_centre(self):
"""Whether we are moving away from [0,0] on the camera.
If we have recorded more than two steps, this property will be
`True` if the most recent point in the history is farther away from
`[0,0]` than the second most recent point. If we have recorded fewer
than 2 points, this property returns None.
"""
if len(self.image_positions) < 2:
return None
else:
return norm(self.image_positions[-1,:]) > norm(self.image_positions[-2])
def move_until_motion_detected(tracker, move, displacement, threshold=10, multipliers=2**np.arange(16), detect_cumulative_motion=False):
"""Move the stage until we can detect motion in the camera.
We move the stage in the direction given by ``displacement`` until the
image has shifted by at least ``threshold`` pixels. The steps will be
given by ``multipliers``, i.e. each time we move to
``displacement * multipliers[i]`` relative to the starting position.
NB we expect that the ``tracker`` object has already been initialised
with ``acquire_template``.
``detect_cumulative_motion`` will use the first point in the tracker as
the point to detect displacement relative to, rather than the last point.
The displacements are always made relative to the last point in the tracker
as it is passed in (i.e. the stage is always moved relative to where it
currently is, but motion detection may be done relative to where the tracker
was initialised). This only matters if the tracker has more than one point
in its history.
The return value `i, m` is the number of moves made, and the largest
multiplier value that was used, i.e. we moved by a total of
`displacement * m`.
"""
displacement = np.array(displacement)
starting_image_position = tracker.image_positions[0 if detect_cumulative_motion else -1, :]
starting_stage_position = tracker.stage_positions[-1, :]
for i, m in enumerate(multipliers):
move(starting_stage_position + displacement * m)
tracker.append_point()
if norm(tracker.image_positions[-1, :] - starting_image_position) >= threshold:
return i + 1, m
raise Exception("Moved the stage by {} but saw no motion.".format(multipliers[-1] * displacement))
def concatenate_tracker_histories(histories):
"""Combine a number of separate tracker history entries into one
A "tracker history" refers to the output of `Tracker.history`, i.e.
it is a tuple of `(stage_positions, image_positions)` with the two
components being a Nx3 and Nx2 `numpy.ndarray` objects respectively.
Given an array of such tuples, we will concatenate the components,
returning a single "tracker history" with the segments concatenated.
Returns: backlash, pixels_per_step, fractional_error
The return value is a tuple of 3 numbers; the estimated backlash (in
motor steps), the ratio of image_position changes to stage_position
(in units of pixels/steps), and an estimate of goodness of fit.
"""
components = zip(*histories)
return tuple(np.concatenate(c, axis=1) for c in components)

94
openflexure_microscope/api/default_extensions/camera_stage_mapping/extension.py Normal file
View file

@ -0,0 +1,94 @@
"""
API extension for stage calibration
This file contains the HTTP API for camera/stage calibration.
"""
from labthings.server.view import View
from labthings.server.find import find_component
from labthings.server.extensions import BaseExtension
from labthings.server.decorators import marshal_task, ThingAction, use_args
from labthings.server import fields
from labthings.core.tasks import taskify
from labthings.core.utilities import get_by_path, set_by_path, create_from_path
from flask import abort
import logging
import time
import numpy as np
import PIL
import io
from .camera_stage_calibration_1d import calibrate_backlash_1d, image_to_stage_displacement_from_1d
from .camera_stage_tracker import Tracker
def update_extension_settings(microscope, settings):
"""Update the stored extension settings dictionary"""
keys = ["extensions","org.openflexure.camera_stage_mapping"]
dictionary = create_from_path(keys)
set_by_path(dictionary, keys, settings)
microscope.update_settings(dictionary)
microscope.save_settings()
def get_extension_settings(microscope):
"""Retrieve the settings for this extension"""
keys = ["extensions","org.openflexure.camera_stage_mapping"]
return get_by_path(microscope.read_settings, keys)
def camera_stage_functions(microscope):
"""Return functions that allow us to interface with the microscope"""
microscope.camera.start_worker() # ensure the worker thread is running, so there is an MJPEG stream
def grab_image():
jpeg = microscope.camera.get_frame()
return np.array(PIL.Image.open(io.BytesIO(jpeg)))
def get_position():
return microscope.stage.position
move = microscope.stage.move_abs
return grab_image, get_position, move
def calibrate_1d(microscope, direction):
"""Move a microscope's stage in 1D, and figure out the relationship with the camera"""
grab_image, get_position, move = camera_stage_functions(microscope)
def wait():
time.sleep(0.2)
tracker = Tracker(grab_image, get_position, settle=wait)
return calibrate_backlash_1d(tracker, move, direction)
@ThingAction
class Calibrate1DView(View):
@use_args({
"direction": fields.List(fields.Float(), required=True, example=[1,0,0])
})
@marshal_task
def post(self, args):
"""Calibrate one axis of the microscope stage against the camera."""
microscope = find_component("org.openflexure.microscope")
direction = np.array(args.get("direction"))
task = taskify(calibrate_1d)(microscope, direction)
return task
csm_extension = BaseExtension("org.openflexure.camera_stage_mapping", version="0.0.1")
csm_extension.add_method(calibrate_1d, "calibrate_1d")
#csm_extension.add_method(calibrate_xy, "calibrate_xy")
csm_extension.add_view(Calibrate1DView, "/calibrate_1d")

237
openflexure_microscope/api/default_extensions/camera_stage_mapping/image_with_location.py Normal file
View file

@ -0,0 +1,237 @@
"""
Image With Location
===================
This datatype supports the various operations that rely on linking a camera to a microscope stage. It is an image
along with the metadata required to relate positions in the image to positions in real life.
To create an `ImageWithLocation`, first put the image data into an `ArrayWithAttrs` and then specify the required
metadata in the `attrs` dictionary. The `pixel_to_sample_matrix` is the only required piece of metadata - the
`datum_pixel` is optional (and if missing, will be assumed to be the central pixel).
A note on coordinate systems
----------------------------
I've tried to stick to two coordinate systems: that used by the stage, generally called a "location", and pixels in an
image.
Images have a "datum pixel", specified in metadata or assumed to be the centre (i.e. pixel (N-1)/2 for a width of N).
This need not be an integer pixel position, but is specified in pixels relative to the [0,0] pixel. When considering
something within an image, *the coordinate system is always relative to pixel [0,0]*, not relative to the datum pixel.
Similarly, the transformation matrix that moves between pixel and stage coordinates uses [0,0] as its origin, not the
datum pixel. However, when considering the displacement between two images, this is usually with respect to the datum
pixels of the images - though we should generally specify this.
We transform between pixel and location coordinate systems with a matrix, the `pixel_to_sample_matrix`. Usually it
is called ``M`` in mathematical expressions. To convert a pixel coordinate to a location, we post-multiply the pixel
coordinate by the matrix, i.e. ``l = p.M`` and to convert the other way we use the inverse of ``M`` so ``p = l.M``
where the dot denotes matrix multiplication using `numpy.dot`.
Note that the calibration matrix is a 4x4 matrix, and the vectors should be ``(x, y, z, 1)`` so that we encode
the absolute position in that matrix, along with scaling and rotation. It's an entirely sensible thing to include
the stage coordinates in metadata as well as the matrix, but it is not needed.
"""
from __future__ import division
from builtins import range
from past.utils import old_div
import numpy as np
from array_with_attrs import ArrayWithAttrs, ensure_attrs
import cv2
#import cv2.cv
from scipy import ndimage
class ImageWithLocation(ArrayWithAttrs):
"""An image, as a numpy array, with attributes to provide location information
This is a functioning `numpy.ndarray` which can store the image in uncompressed format.
We require that the `attrs` dictionary (defined by `ArrayWithAttrs`) contains keys
that we use to store the crucial mapping from pixels in the image to position in the
sample.
"""
# def __array_finalize__(self, obj):
# """Ensure that the object is a properly set-up ImageWithLocation"""
# ArrayWithAttrs.__array_finalize__(self, obj) # Ensure we have self.attrs
def __getitem__(self, item):
"""Update the metadata when we extract a slice"""
try:
# Handle specially the case where we are extracting a 2D region of the image, i.e. the first and second
# indices are slices. We test for that here - and do it in a try: except block so that if, for example,
# item is not indexable,
assert isinstance(item[0], slice), "First index was not a slice"
assert isinstance(item[1], slice), "Second index was not a slice"
start = np.array([item[i].start for i in range(2)])
start = np.where(start == np.array(None), 0, start) # missing start points are equivalent to zero
step = np.array([item[i].step for i in range(2)])
step = np.where(step == np.array(None), 1, step) # missing step is equivalent to step==1
except:
# If the above doesn't work, assume we're not dealing with a 2D slice and give up.
return super(ImageWithLocation, self).__getitem__(item) # pass it on up
out = super(ImageWithLocation, self).__getitem__(item) # retrieve the slice
out.datum_pixel -= start # adjust the datum pixel so it refers to the same part of the image
# Next, we adjust the constant part of the pixel-sample matrix so pixels stay in the same place
location_shift = np.dot(ensure_3d(start), self.pixel_to_sample_matrix[:3,:3])
out.pixel_to_sample_matrix[3,:3] += location_shift
if not np.all(step == 1):
# if we're downsampling, remember to scale datum_pixel accordingly
out.datum_pixel = old_div(out.datum_pixel, step)
# Scale the pixel-to-sample matrix if we've got a non-unity step in the slice
# I don't understand why I can't do this with slicing, but it all goes wrong...
for i in range(2):
out.pixel_to_sample_matrix[i, :3] *= step[i]
return out
def pixel_to_location(self, pixel):
"""Return the location in the sample of the given pixel.
NB this returns a 3D location, including Z."""
p = ensure_2d(pixel)
l = np.dot(np.array([p[0], p[1], 0, 1]), self.pixel_to_sample_matrix)
return l[:3]
def location_to_pixel(self, location, check_bounds=False, z_tolerance=np.infty):
"""Return the pixel coordinates of a given location in the sample.
location : numpy.ndarray
A 2- or 3- element numpy array representing sample position, in units of distance.
check_bounds : bool, optional (default False)
If this is True, raise an exception if the pixel is not in the image.
z_tolerance : float, optional (defaults to infinity)
If we are checking the bounds, make sure the sample location is within this distance of the image's Z
position. The default is to allow any distance.
Returns : numpy.ndarray
A 2- or 3- element position, to match the size of location passed in.
"""
l = ensure_2d(location)
l = l[:2]-self.pixel_to_sample_matrix[3,:2]
p = np.dot(l, np.linalg.inv(self.pixel_to_sample_matrix[:2,:2]))
if check_bounds:
assert np.all(0 <= p[0:2]), "The location was not within the image"
assert np.all(p[0:2] <= self.shape[0:2]), "The location was not within the image"
assert np.abs(p[2]) < z_tolerance, "The location was too far away from the plane of the image"
if len(location) == 2:
return p[:2]
else:
return p[:3]
def feature_at(self, centre_position, size=(100,100), set_datum_to_centre=True):
"""Return a thumbnail cropped out of this image, centred on a particular pixel position.
This is simply a convenience method that saves typing over the usual slice syntax. Below are two equivalent
ways of extracting a thumbnail:
pos = (240,320)
size = (100,100)
thumbnail = image[pos[0] - size[0]/2:pos[0] + size[0]/2, pos[1] - size[1]/2:pos[1] + size[1]/2, ...]
thumbnail2 = image.feature_at(pos, size)
thumbnail3 = image[190:290 270:370]
``centre_position`` and ``size`` should be two-element tuples, but the intention is that this code will cope
gracefully with floating-point values.
NB the datum pixel of the returned image will be set to its centre, not the datum position of the original image
by default. Give the argument ``set_datum_to_centre=False`` to disable this behaviour.
"""
try:
float(centre_position[0])
float(centre_position[1])
float(size[0])
float(size[1])
except:
raise IndexError("Error: arguments of feature_at were invalid: {}, {}".format(centre_position, size))
pos = centre_position
# For now, rely on numpy to complain if the feature is outside the image. May do bound-checking at some point.
# If so, we might need to think carefully about the datum pixel of the resulting image.
thumb = self[pos[0] - old_div(size[0],2):pos[0] + old_div(size[0],2), pos[1] - old_div(size[1],2):pos[1] + old_div(size[1],2), ...]
if set_datum_to_centre:
thumb.datum_pixel = (old_div(size[0],2), old_div(size[1],2)) # Make the datum point of the new image its centre.
return thumb
def downsample(self, n):
"""Return a view of the image, downsampled (sliced with a non-unity step).
In the future, an optional argument to this function may take means of blocks of the images to improve signal
to noise. Currently it just decimates (i.e. throws away rows and columns).
"""
assert n > 0, "The downsampling factor must be an integer greater than 0"
return self[::int(n), ::int(n), ...] # The slicing code handles updating metadata
@property
def datum_pixel(self):
"""The pixel that nominally corresponds to where the image "is".
Usually the datum pixel is the central pixel, and if the metadata required is not present,
we will silently assume that this is the case.
"""
datum = self.attrs.get('datum_pixel', old_div((np.array(self.shape[:2]) - 1),2))
assert len(datum) == 2, "The datum pixel didn't have length 2!"
return datum
@datum_pixel.setter
def datum_pixel(self, datum):
assert len(datum) == 2, "The datum pixel didn't have length 2!"
self.attrs['datum_pixel'] = datum
@property
def datum_location(self):
"""The location in the sample of the datum pixel"""
return self.pixel_to_location(self.datum_pixel)
@property
def pixel_to_sample_matrix(self):
"""The matrix that maps from pixel coordinates to sample coordinates.
np.dot(p, M) yields a location for the given pixel, where p is [x,y,0,1] and M is this matrix. The location
given will be 4 elements long, and will have 1 as the final element.
"""
M = self.attrs['pixel_to_sample_matrix']
assert M.shape == (4, 4), "The pixel-to-sample matrix is the wrong shape!"
assert M.dtype.kind == "f", "The pixel-to-sample matrix is not floating point!"
return M
@pixel_to_sample_matrix.setter
def pixel_to_sample_matrix(self, M):
M = np.asanyarray(M) #ensure it's an ndarray subclass
assert M.shape == (4, 4), "The pixel-to-sample matrix must be 4x4!"
assert M.dtype.kind == "f", "The pixel-to-sample matrix must be floating point!"
self.attrs['pixel_to_sample_matrix'] = M
#TODO: split the data type out of this module and put it somewhere sensible
def add_location_metadata(image, pixel_to_sample_matrix, datum_pixel=None):
"""Wrap an image if needed, and set its pixel to sample matrix."""
awa = ensure_attrs(image) # if needed, convert the image to an ArrayWithAttrs
awa.attrs['pixel_to_sample_matrix'] = pixel_to_sample_matrix
if datum_pixel is not None:
awa.attrs['datum_pixel'] = datum_pixel
return awa
def datum_pixel(image):
"""Get the datum pixel of an image - if no property is present, assume the central pixel."""
try:
return np.array(image.datum_pixel)
except:
return (np.array(image.shape[:2]) - 1) / 2.
def ensure_3d(vector):
"""Make sure a vector has 3 elements, appending a zero if needed."""
if len(vector) == 3:
return np.array(vector)
elif len(vector) == 2:
return np.array([vector[0], vector[1], 0])
else:
raise ValueError("Tried to ensure a vector was 3D, but it had neither 2 nor 3 elements!")
def ensure_2d(vector):
"""Make sure a vector has 3 elements, appending a zero if needed."""
if len(vector) == 2:
return np.array(vector)
elif len(vector) == 3:
return np.array(vector[:2])
else:
raise ValueError("Tried to ensure a vector was 2D, but it had neither 2 nor 3 elements!")