Source code for shampoo.focus

from __future__ import (absolute_import, division, print_function,

from .reconstruction import ReconstructedWave, unwrap_phase

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN

__all__ = ['cluster_focus_peaks', 'find_focus_plane', 'locate_specimens']

[docs]def cluster_focus_peaks(xyz, eps=5, min_samples=3): """ Use DBSCAN to identify single particles through multiple focus planes. Parameters ---------- xyz : `~numpy.ndarray` Matrix of (x, y, z) positions for each peak detected eps : float Passed to the ``eps`` argument of `~sklearn.cluster.DBSCAN` min_samples : int Passed to the ``min_samples`` argument of `~sklearn.cluster.DBSCAN` Returns ------- labels : `~numpy.ndarray` List of cluster labels for each peak. Labels of `-1` signify noise points. """ positions = xyz.copy() # Compress distances in the z-axis positions[:, 2] /= 10 db = DBSCAN(eps=eps, min_samples=min_samples).fit(positions) labels = db.labels_ return labels
[docs]def find_focus_plane(roi_cube, focus_on='amplitude', plot=False): """ Find focus plane in a cube of reconstructed waves at different propagation distances. Uses the autofocus method of Dubois et al. 2006 who showed that the integral over the image plane of the amplitude of the reconstructed wave is minimum at the focal plane for a pure amplitude object [1]_. This will only work for small cubes centered on the specimen. .. [1] Parameters ---------- roi_cube : `~numpy.ndarray` Reconstructed waves at ``N`` propagation distances with ``M`` by ``M`` pixels, with a shape of ``(N, M, M)`` focus_on : {"amplitude", "phase"} (optional) Focus on the phase or amplitude? plot : bool (optional) Make plots of the integral of the amplitude of the reconstructed wave as a function of distance. Default is False. Returns ------- focus_index : int Index of the z-plane that is in focus significance : float Significance of a detection of a particle in the ROI cube, measured by the significance of the negative peak in the derivative of the sum of the unwrapped phase in the image plane, with respect to propagation distance. For example, if ``significance < 3`` the detection of a specimen may be legitimate. """ if focus_on == 'amplitude': extremum = np.argmin elif focus_on == 'phase': extremum = np.argmax else: raise ValueError('The `focus_on` kwarg must be either "phase" or ' '"amplitude".') # Following Equation 9, 10 of Dubois et al. 2006: integral_abs_wave = np.sum(np.abs(roi_cube), axis=(1, 2)) focus_index = extremum(integral_abs_wave) # Do a similar integral on the unwrapped phase. The phase changes # most rapidly on a source near focus, so the derivative wrt propagation # distance of the phase integrated in space has a *minimum* near focus integral_phase_wave = np.sum([unwrap_phase(roi_cube[i, ...]) for i in range(roi_cube.shape[0])], axis=(1, 2)) d_int_phase = np.diff(integral_phase_wave) # Measure significance of detected focus by taking the median normalized, # standard deviation normalized derivative of phase with respect to # propagation distance. This quantity becomes significantly negative for # real specimens. For example, significance < -3 may yield real particles. significance = np.min(d_int_phase / np.median(d_int_phase) / d_int_phase.std()) if plot: plt.figure() plt.plot(range(roi_cube.shape[0]), (integral_abs_wave - integral_abs_wave.mean())/integral_abs_wave.std(), label='intensity') plt.plot(range(roi_cube.shape[0]-1), (d_int_phase - d_int_phase.mean())/d_int_phase.std(), label='d(phase)/d(prop dist)') plt.xlabel('Propagation distance index') plt.title('Significance: {0:.4f}'.format(significance)) plt.axvline(focus_index, ls='--') plt.legend() return focus_index, significance
def _correct_limits(minimum, maximum, axis_range, edge): if minimum < axis_range: minimum = axis_range if maximum > edge - axis_range: maximum = edge - axis_range return minimum, maximum, axis_range
[docs]def locate_specimens(wave_cube, positions, labels, distances, plots=False): """ Identify the (x, y, z) coordinates of a specimen. Parameters ---------- wave_cube : `~numpy.ndarray` Cube of reconstructed waves positions : `~numpy.ndarray` (x,y,z) positions of objects detected by the blob finder labels : `~numpy.ndarray` Clustering labels for each (x,y,z) coordinate, identifying groups of positions, i.e., single particles detected at multiple z-planes distances : `~numpy.ndarray` Propagation distances, same length as the first axis of ``complex_cube`` Returns ------- specimen_coordinates : `~numpy.ndarray` (x, y, z) coordinates of each detected specimen specimen_coordinates : `~numpy.ndarray` Significance of each specimen detection. See docs of `~shampoo.focus.find_focus_plane` for hints on how to interpret the significance quantity. """ specimen_coordinates = [] specimen_significance = [] for l in set(labels): n_points = np.count_nonzero(labels == l) if l != -1 and n_points > 3: xmedian = np.median(positions[labels == l, 0]) ymedian = np.median(positions[labels == l, 1]) xmin, ymin, zmin_d = np.min(positions[labels == l, :], axis=0) xmax, ymax, zmax_d = np.max(positions[labels == l, :], axis=0) zmin = np.argmin(np.abs(zmin_d - distances)) zmax = np.argmin(np.abs(zmax_d - distances)) x_range = y_range = 2 z_range = zmax - zmin xmin, xmax, x_range = _correct_limits(xmin, xmax, x_range, wave_cube.shape[1]) ymin, ymay, y_range = _correct_limits(ymin, ymax, y_range, wave_cube.shape[2]) zmin, zmaz, z_range = _correct_limits(zmin, zmax, z_range, wave_cube.shape[0]) # Make a reconstructed wave cube centered on the region of interest roi_cube = wave_cube[zmin - z_range:zmax + z_range, xmin - x_range:xmax + x_range, ymin - y_range:ymax + y_range] # Using this cropped cube centered on the ROI, find the best focus focus_ind_minus_margin, significance = find_focus_plane(roi_cube, plot=plots) focus_ind = focus_ind_minus_margin + zmin - z_range specimen_coordinates.append([xmedian, ymedian, distances[focus_ind]]) specimen_significance.append(significance) if plots: focused_wave = ReconstructedWave(wave_cube[focus_ind, ...]) fig, ax = focused_wave.plot(phase=True) thetas = np.linspace(0, 2*np.pi, 30) r = 20 ax.plot(r*np.cos(thetas) + ymedian, r*np.sin(thetas) + xmedian, lw=3, color='r') return np.array(specimen_coordinates), np.array(specimen_significance)