diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 8356d0be..e42ced66 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -8,7 +8,11 @@ version: 2 build: os: ubuntu-22.04 tools: - python: "3.10" + python: + version: 3.10 + install: + - method: pip + - path: . # Build documentation in the "docs/" directory with Sphinx sphinx: diff --git a/changelog.rst b/changelog.rst index b603bd1d..657641e7 100644 --- a/changelog.rst +++ b/changelog.rst @@ -3,10 +3,16 @@ Changelog Last change: 22-FEB-2024 MTS -0.6.9 ------ +0.6.9 - 0.6.10 +-------------- - Added the option to draw polygon picks in Picasso: Render - Save pick properties in Picasso: Render saves areas of picked regions in nm^2 +- Calibration .yaml file saves number of frames and step size in nm +- ``picasso.lib.merge_locs`` function can merge localizations from multiple files +- Mask dialog in Picasso: Render saves .png mask files +- Mask dialog in Picasso: Render allows to save .png with the blurred image +- Picasso: Localize - added the option to save the current view as a .png file +- Picasso:Render - functions related to picking moved to ``picasso.lib`` and ``picasso.postprocess`` 0.6.6 - 0.6.8 ------------- diff --git a/picasso/__main__.py b/picasso/__main__.py index d48e8533..d14bef8e 100644 --- a/picasso/__main__.py +++ b/picasso/__main__.py @@ -165,15 +165,17 @@ def _hdf2csv(path): paths = glob(path) if paths: import os.path + from .io import load_locs for path in _tqdm(paths): base, ext = os.path.splitext(path) if ext == ".hdf5": print("Converting {}".format(path)) out_path = base + ".csv" - locs = pd.read_hdf(path) + locs = load_locs(path)[0] + df = pd.DataFrame(locs) print("A total of {} rows loaded".format(len(locs))) - locs.to_csv(out_path, sep=",", encoding="utf-8") + df.to_csv(out_path, sep=",", encoding="utf-8") print("Complete.") diff --git a/picasso/gui/localize.py b/picasso/gui/localize.py index acc5de25..6a3df19a 100644 --- a/picasso/gui/localize.py +++ b/picasso/gui/localize.py @@ -1804,11 +1804,12 @@ def export_current(self): self, "Save image", out_path, filter="*.png;;*.tif" ) if path: - qimage = QtGui.QImage(self.view.viewport().size(), QtGui.QImage.Format_RGB32) + qimage = QtGui.QImage(self.view.size(), QtGui.QImage.Format_RGB32) painter = QtGui.QPainter(qimage) self.view.render(painter) painter.end() qimage.save(path) + # self.view.scene().save_image(path) self.view.setMinimumSize(1, 1) def save_locs(self, path): diff --git a/picasso/gui/render.py b/picasso/gui/render.py index 18f93323..d2419a11 100644 --- a/picasso/gui/render.py +++ b/picasso/gui/render.py @@ -3582,7 +3582,7 @@ def __init__(self, window): self.save_mask_button.clicked.connect(self.save_mask) mask_grid.addWidget(self.save_mask_button, 5, 1) - self.save_blur_button = QtWidgets.QPushButton("Save Blur") + self.save_blur_button = QtWidgets.QPushButton("Save Blurred") self.save_blur_button.setEnabled(False) self.save_blur_button.setFocusPolicy(QtCore.Qt.NoFocus) self.save_blur_button.clicked.connect(self.save_blur) @@ -3661,8 +3661,8 @@ def save_mask(self): ) if path: np.save(path, self.mask) - png_path = base + "_mask" + ".png" - plt.imsave(png_path, self.mask, cmap='gray') + png_path = path.replace(".npy", ".png") + plt.imsave(png_path, self.mask, cmap="gray") def save_blur(self): """Save blurred image to a png. format.""" @@ -5386,11 +5386,7 @@ class View(QtWidgets.QLabel): Finds group color index for each localization get_index_blocks(channel) Calls self.index_locs if not calculated earlier - get_pick_polygon_corners(pick) - Returns X and Y coordinates of a pick polygon - get_pick_rectangle_corners(start_x, start_y, end_x, end_y, width) - Finds the positions of a rectangular pick's corners - get_pick_rectangle_polygon(start_x, start_y, end_x, end_y, width) + get_pick_polygon(start_x, start_y, end_x, end_y, width) Finds a PyQt5 object used for drawing a rectangular pick get_render_kwargs() Returns a dictionary to be used for the kwargs of render.render @@ -5432,7 +5428,7 @@ class View(QtWidgets.QLabel): pan_relative(dy, dx) Moves viewport by a given relative distance pick_areas() - Finds the areas of all current picks in nm^2. + Finds the areas of all current picks in um^2. pick_message_box(params) Returns a message box for selecting picks pick_similar() @@ -6050,7 +6046,7 @@ def dbscan(self): "Input Dialog", "Enter suffix", QtWidgets.QLineEdit.Normal, - "_clustered", + "_dbscan", ) if ok: for channel in range(len(self.locs_paths)): @@ -6063,7 +6059,7 @@ def dbscan(self): path, ext = QtWidgets.QFileDialog.getSaveFileName( self, "Save clustered locs", - self.locs_paths[channel].replace(".hdf5", "_clustered.hdf5"), + self.locs_paths[channel].replace(".hdf5", "_dbscan.hdf5"), filter="*.hdf5", ) if path: @@ -6118,7 +6114,7 @@ def _dbscan(self, channel, path, params): status.close() if save_centers: status = lib.StatusDialog("Calculating cluster centers", self) - path = path.replace(".hdf5", "_cluster_centers.hdf5") + path = path.replace(".hdf5", "_centers.hdf5") centers = clusterer.find_cluster_centers(locs, pixelsize=pixelsize) io.save_locs(path, centers, self.infos[channel] + [dbscan_info]) status.close() @@ -6142,7 +6138,7 @@ def hdbscan(self): "Input Dialog", "Enter suffix", QtWidgets.QLineEdit.Normal, - "_clustered", + "_hdbscan", ) if ok: for channel in range(len(self.locs_paths)): @@ -6157,7 +6153,7 @@ def hdbscan(self): "Save clustered locs", self.locs_paths[channel].replace( ".hdf5", - "_clustered.hdf5" + "_hdbscan.hdf5" ), filter="*.hdf5", ) @@ -6215,7 +6211,7 @@ def _hdbscan(self, channel, path, params): status.close() if save_centers: status = lib.StatusDialog("Calculating cluster centers", self) - path = path.replace(".hdf5", "_cluster_centers.hdf5") + path = path.replace(".hdf5", "_centers.hdf5") centers = clusterer.find_cluster_centers(locs, pixelsize=pixelsize) io.save_locs(path, centers, self.infos[channel] + [hdbscan_info]) status.close() @@ -6359,7 +6355,7 @@ def _smlm_clusterer(self, channel, path, params): # save cluster centers if params[-3]: status = lib.StatusDialog("Calculating cluster centers", self) - path = path.replace(".hdf5", "_cluster_centers.hdf5") + path = path.replace(".hdf5", "_centers.hdf5") centers = clusterer.find_cluster_centers(clustered_locs, pixelsize) io.save_locs(path, centers, info) status.close() @@ -6460,50 +6456,7 @@ def dragEnterEvent(self, event): else: event.ignore() - def get_pick_polygon_corners(self, pick): - """Returns X and Y coordinates of a pick polygon. - - Returns None, None if the pick is not a closed polygon.""" - - if len(pick) < 3 or pick[0] != pick[-1]: - return None, None - else: - X = [_[0] for _ in pick] - Y = [_[1] for _ in pick] - return X, Y - - def get_pick_rectangle_corners( - self, start_x, start_y, end_x, end_y, width - ): - """ - Finds the positions of corners of a rectangular pick. - Rectangular pick is defined by: - [(start_x, start_y), (end_x, end_y)] - and its width. (all values in pixels) - - Returns - ------- - tuple - Contains corners' x and y coordinates in two lists - """ - - if end_x == start_x: - alpha = np.pi / 2 - else: - alpha = np.arctan((end_y - start_y) / (end_x - start_x)) - dx = width * np.sin(alpha) / 2 - dy = width * np.cos(alpha) / 2 - x1 = float(start_x - dx) - x2 = float(start_x + dx) - x4 = float(end_x - dx) - x3 = float(end_x + dx) - y1 = float(start_y + dy) - y2 = float(start_y - dy) - y4 = float(end_y + dy) - y3 = float(end_y - dy) - return [x1, x2, x3, x4], [y1, y2, y3, y4] - - def get_pick_rectangle_polygon( + def get_pick_polygon( self, start_x, start_y, end_x, end_y, width, return_most_right=False ): """ @@ -6515,7 +6468,7 @@ def get_pick_rectangle_polygon( QtGui.QPolygonF """ - X, Y = self.get_pick_rectangle_corners( + X, Y = lib.get_pick_rectangle_corners( start_x, start_y, end_x, end_y, width ) p = QtGui.QPolygonF() @@ -6618,7 +6571,7 @@ def draw_picks(self, image): painter.drawLine(start_x, start_y, end_x, end_y) # draw a rectangle - polygon, most_right = self.get_pick_rectangle_polygon( + polygon, most_right = self.get_pick_polygon( start_x, start_y, end_x, end_y, w, return_most_right=True ) painter.drawPolygon(polygon) @@ -6695,7 +6648,7 @@ def draw_rectangle_pick_ongoing(self, image): # convert from camera units to display units w *= self.width() / self.viewport_width() - polygon = self.get_pick_rectangle_polygon( + polygon = self.get_pick_polygon( self.rectangle_pick_start_x, self.rectangle_pick_start_y, self.rectangle_pick_current_x, @@ -7070,13 +7023,13 @@ def move_to_pick(self): xc = np.mean([xs, xe]) yc = np.mean([ys, ye]) w = self.window.tools_settings_dialog.pick_width.value() - X, Y = self.get_pick_rectangle_corners(xs, ys, xe, ye, w) + X, Y = lib.get_pick_rectangle_corners(xs, ys, xe, ye, w) x_min = min(X) - (0.2 * (xc - min(X))) x_max = max(X) + (0.2 * (max(X) - xc)) y_min = min(Y) - (0.2 * (yc - min(Y))) y_max = max(Y) + (0.2 * (max(Y) - yc)) elif self._pick_shape == "Polygon": - X, Y = self.get_pick_polygon_corners(self._picks[pick_no]) + X, Y = lib.get_pick_polygon_corners(self._picks[pick_no]) x_min = min(X) - 0.2 * (max(X) - min(X)) x_max = max(X) + 0.2 * (max(X) - min(X)) y_min = min(Y) - 0.2 * (max(Y) - min(Y)) @@ -8573,7 +8526,7 @@ def get_index_blocks(self, channel, fast_render=False): @check_pick def pick_areas(self): - """Finds the areas of all current picks in nm^2. + """Finds the areas of all current picks in um^2. Returns ------- @@ -8584,25 +8537,15 @@ def pick_areas(self): if self._pick_shape == "Circle": d = self.window.tools_settings_dialog.pick_diameter.value() r = d / 2 - areas = np.ones(len(self._picks)) * np.pi * r ** 2 + areas = lib.pick_areas_circle(self._picks, r) elif self._pick_shape == "Rectangle": w = self.window.tools_settings_dialog.pick_width.value() - areas = np.zeros(len(self._picks)) - for i, pick in enumerate(self._picks): - (xs, ys), (xe, ye) = pick - areas[i] = w * np.sqrt((xe - xs) ** 2 + (ye - ys) ** 2) + areas = lib.pick_areas_rectangle(self._picks, w) elif self._pick_shape == "Polygon": - areas = np.zeros(len(self._picks)) - for i, pick in enumerate(self._picks): - if len(pick) < 3 or pick[0] != pick[-1]: # not a closed polygon - areas[i] = 0 - continue - X, Y = self.get_pick_polygon_corners(pick) - areas[i] = lib.polygon_area(X, Y) - areas = areas[areas > 0] # remove open polygons - + areas = lib.pick_areas_polygon(self._picks) + pixelsize = self.window.display_settings_dlg.pixelsize.value() - areas *= pixelsize ** 2 + areas *= (pixelsize * 1e-3) ** 2 # convert to um^2 return areas @check_picks @@ -8698,112 +8641,52 @@ def picked_locs( Parameters ---------- channel : int - Channel of locs to be processed + Channel of locs to be processed. add_group : boolean (default=True) True if group id should be added to locs. Each pick will be - assigned a different id + assigned a different id. fast_render : boolean If True, takes self.locs, i.e. after randomly sampling a - fraction of self.all_locs. If False, takes self.all_locs + fraction of self.all_locs. If False, takes self.all_locs. Returns ------- - list - List of np.recarrays, each containing locs from one pick + picked_locs : list of np.recarrays + List of np.recarrays, each containing locs from one pick. """ if len(self._picks): - picked_locs = [] + # initialize progress dialog progress = lib.ProgressDialog( "Creating localization list", 0, len(self._picks), self ) progress.set_value(0) + + # extract localizations to pick from + if fast_render: + locs = self.locs[channel].copy() + else: + locs = self.all_locs[channel].copy() + + # find pick size (radius or width) if self._pick_shape == "Circle": d = self.window.tools_settings_dialog.pick_diameter.value() - r = d / 2 - index_blocks = self.get_index_blocks( - channel, fast_render=fast_render - ) - for i, pick in enumerate(self._picks): - x, y = pick - block_locs = postprocess.get_block_locs_at( - x, y, index_blocks - ) - group_locs = lib.locs_at(x, y, block_locs, r) - if add_group: - group = i * np.ones(len(group_locs), dtype=np.int32) - group_locs = lib.append_to_rec( - group_locs, group, "group" - ) - group_locs.sort(kind="mergesort", order="frame") - picked_locs.append(group_locs) - progress.set_value(i + 1) + pick_size = d / 2 elif self._pick_shape == "Rectangle": - w = self.window.tools_settings_dialog.pick_width.value() - if fast_render: - channel_locs = self.locs[channel] - else: - channel_locs = self.all_locs[channel] - for i, pick in enumerate(self._picks): - (xs, ys), (xe, ye) = pick - X, Y = self.get_pick_rectangle_corners(xs, ys, xe, ye, w) - x_min = min(X) - x_max = max(X) - y_min = min(Y) - y_max = max(Y) - group_locs = channel_locs[channel_locs.x > x_min] - group_locs = group_locs[group_locs.x < x_max] - group_locs = group_locs[group_locs.y > y_min] - group_locs = group_locs[group_locs.y < y_max] - group_locs = lib.locs_in_rectangle(group_locs, X, Y) - # store rotated coordinates in x_rot and y_rot - angle = 0.5 * np.pi - np.arctan2((ye - ys), (xe - xs)) - x_shifted = group_locs.x - xs - y_shifted = group_locs.y - ys - x_pick_rot = x_shifted * np.cos( - angle - ) - y_shifted * np.sin(angle) - y_pick_rot = x_shifted * np.sin( - angle - ) + y_shifted * np.cos(angle) - group_locs = lib.append_to_rec( - group_locs, x_pick_rot, "x_pick_rot" - ) - group_locs = lib.append_to_rec( - group_locs, y_pick_rot, "y_pick_rot" - ) - if add_group: - group = i * np.ones(len(group_locs), dtype=np.int32) - group_locs = lib.append_to_rec( - group_locs, group, "group" - ) - group_locs.sort(kind="mergesort", order="frame") - picked_locs.append(group_locs) - progress.set_value(i + 1) - elif self._pick_shape == "Polygon": - if fast_render: - channel_locs = self.locs[channel] - else: - channel_locs = self.all_locs[channel] - for i, pick in enumerate(self._picks): - X, Y = self.get_pick_polygon_corners(pick) - if X is None: - progress.set_value(i + 1) - continue - group_locs = channel_locs[channel_locs.x > min(X)] - group_locs = group_locs[group_locs.x < max(X)] - group_locs = group_locs[group_locs.y > min(Y)] - group_locs = group_locs[group_locs.y < max(Y)] - group_locs = lib.locs_in_polygon(group_locs, X, Y) - if add_group: - group = i * np.ones(len(group_locs), dtype=np.int32) - group_locs = lib.append_to_rec( - group_locs, group, "group" - ) - group_locs.sort(kind="mergesort", order="frame") - picked_locs.append(group_locs) - progress.set_value(i + 1) - + pick_size = self.window.tools_settings_dialog.pick_width.value() + else: + pick_size = None + + # pick localizations + picked_locs = postprocess.picked_locs( + locs, + self.infos[channel], + self._picks, + self._pick_shape, + pick_size=pick_size, + add_group=add_group, + callback=progress.set_value, + ) return picked_locs def remove_picks(self, position): @@ -8832,7 +8715,7 @@ def remove_picks(self, position): y = np.array([y]) for pick in self._picks: (start_x, start_y), (end_x, end_y) = pick - X, Y = self.get_pick_rectangle_corners( + X, Y = lib.get_pick_rectangle_corners( start_x, start_y, end_x, end_y, width ) # do not check if rectangle has no size @@ -9259,9 +9142,13 @@ def save_picked_locs(self, path, channel): # save picked locs with .yaml if locs is not None: + areas = self.pick_areas() pick_info = { "Generated by": "Picasso Render : Pick", "Pick Shape": self._pick_shape, + "Pick Areas (um^2)": [float(_) for _ in areas], + "Total Picked Area (um^2)": float(np.sum(areas)), + "Number of picks": len(self._picks), } if self._pick_shape == "Circle": d = self.window.tools_settings_dialog.pick_diameter.value() @@ -9296,10 +9183,14 @@ def save_picked_locs_multi(self, path): # save locs = locs.view(np.recarray) if locs is not None: - d = self.window.tools_settings_dialog.pick_diameter.value() + areas = self.pick_areas() pick_info = { "Generated by:": "Picasso Render : Pick", "Pick Shape:": self._pick_shape, + "Pick Areas (um^2)": [float(_) for _ in areas], + "Total Picked Area (um^2)": float(np.sum(areas)), + "Number of picks": len(self._picks), + } if self._pick_shape == "Circle": d = self.window.tools_settings_dialog.pick_diameter.value() @@ -9363,7 +9254,7 @@ def save_pick_properties(self, path, channel): ) # add the area of the picks to the properties areas = self.pick_areas() - pick_props = lib.append_to_rec(pick_props, areas, "pick_area_nm2") + pick_props = lib.append_to_rec(pick_props, areas, "pick_area_um2") progress.close() # QPAINT estimate of number of binding sites n_units = self.window.info_dialog.calculate_n_units(dark) @@ -11195,6 +11086,8 @@ def export_current_info(self, path): "Scalebar Length (nm)": d.scalebar.value(), "Localizations Loaded": self.view.locs_paths, "Colors": colors, + "Display pixel size (nm)": d.disp_px_size.value(), + "Min. blur (cam. px)": d.min_blur_width.value(), } io.save_info(path, [info]) @@ -11215,6 +11108,7 @@ def export_complete(self): viewport = [(0, 0), (movie_height, movie_width)] qimage = self.view.render_scene(cache=False, viewport=viewport) qimage.save(path) + self.export_current_info(path) def export_txt(self): """ diff --git a/picasso/lib.py b/picasso/lib.py index d981c5e4..bae337bc 100644 --- a/picasso/lib.py +++ b/picasso/lib.py @@ -14,6 +14,7 @@ from lmfit import Model as _Model from numpy.lib.recfunctions import append_fields as _append_fields from numpy.lib.recfunctions import drop_fields as _drop_fields +from numpy.lib.recfunctions import stack_arrays as _stack_arrays import collections as _collections import glob as _glob import os.path as _ospath @@ -27,6 +28,8 @@ class ProgressDialog(QtWidgets.QProgressDialog): + """ProgressDialog displays a progress dialog with a progress bar.""" + def __init__(self, description, minimum, maximum, parent): super().__init__( description, @@ -56,6 +59,8 @@ def closeEvent(self, event): class StatusDialog(QtWidgets.QDialog): + """StatusDialog displays the description string in a dialog.""" + def __init__(self, description, parent): super(StatusDialog, self).__init__(parent, QtCore.Qt.CustomizeWindowHint) _dialogs.append(self) @@ -82,6 +87,8 @@ def __init__(self, *args, **kwargs): def cancel_dialogs(): + """Closes all open dialogs.""" + dialogs = [_ for _ in _dialogs] for dialog in dialogs: if isinstance(dialog, ProgressDialog): @@ -99,6 +106,21 @@ def cumulative_exponential(x, a, t, c): def calculate_optimal_bins(data, max_n_bins=None): + """Calculates the optimal bins for display. + + Parameters + ---------- + data : numpy.1darray + Data to be binned. + max_n_bins : int (default=None) + Maximum number of bins. + + Returns + ------- + bins : numpy.1darray + Bins for display. + """ + iqr = _np.subtract(*_np.percentile(data, [75, 25])) bin_size = 2 * iqr * len(data) ** (-1 / 3) if data.dtype.kind in ("u", "i") and bin_size < 1: @@ -111,13 +133,31 @@ def calculate_optimal_bins(data, max_n_bins=None): return None if max_n_bins and n_bins > max_n_bins: n_bins = max_n_bins - return _np.linspace(bin_min, data.max(), n_bins) + bins = _np.linspace(bin_min, data.max(), n_bins) + return bins def append_to_rec(rec_array, data, name): + """Appends a new column to the existing np.recarray. + + Parameters + ---------- + rec_array : np.rec.array + Recarray to which the new column is appended. + data : np.1darray + Data to be appended. + name : str + Name of the new column. + + Returns + ------- + rec_array : np.rec.array + Recarray with the new column. + """ + if hasattr(rec_array, name): rec_array = remove_from_rec(rec_array, name) - return _append_fields( + rec_array = _append_fields( rec_array, name, data, @@ -127,7 +167,53 @@ def append_to_rec(rec_array, data, name): ) return rec_array + +def merge_locs(locs_list, increment_frames=True): + """Merges localization lists into one file. Can increment frames + to avoid overlapping frames. + + Parameters + ---------- + locs_list : list of np.rec.arrays + List of localization lists to be merged. + increment_frames : bool (default=True) + If True, increments frames of each localization list by the + maximum frame number of the previous localization list. Useful + when the localization lists are from different movies but + represent the same stack. + + Returns + locs : np.rec.array + Merged localizations. + """ + + if increment_frames: + last_frame = 0 + for i, locs in enumerate(locs_list): + locs["frame"] += last_frame + last_frame = locs["frame"][-1].max() + locs_list[i] = locs + locs = _stack_arrays(locs_list, usemask=False, asrecarray=True) + return locs + + def ensure_sanity(locs, info): + """Ensures that localizations are within the image dimensions + and have positive localization precisions. + + Parameters + ---------- + locs : np.rec.array + Localizations list. + info : list of dicts + Localization metadata. + + Returns + ------- + locs : np.rec.array + Localizations that pass the sanity checks. + """ + # no inf or nan: locs = locs[ _np.all( @@ -146,40 +232,55 @@ def ensure_sanity(locs, info): def is_loc_at(x, y, locs, r): + """Checks if localizations are at position (x, y) within radius r. + + Parameters + ---------- + x : float + x-coordinate of the position. + y : float + y-coordinate of the position. + locs : np.rec.array + Localizations list. + r : float + Radius. + + Returns + ------- + is_picked : np.ndarray + Boolean array indicating if localization is at position. + """ + dx = locs.x - x dy = locs.y - y r2 = r**2 - return dx**2 + dy**2 < r2 + is_picked = dx**2 + dy**2 < r2 + return is_picked def locs_at(x, y, locs, r): - is_picked = is_loc_at(x, y, locs, r) - return locs[is_picked] - + """Returns localizations at position (x, y) within radius r. -def polygon_area(X, Y): - """Finds the area of a polygon defined by corners X and Y. - Parameters ---------- - X : numpy.1darray - x-coordinates of the polygon corners. - Y : numpy.1darray - y-coordinates of the polygon corners. - - Returns + x : float + x-coordinate of the position. + y : float + y-coordinate of the position. + locs : np.rec.array + Localizations list. + r : float + Radius. + + Returns ------- - area : float - Area of the polygon. + picked_locs : np.rec.array + Localizations at position. """ - n_corners = len(X) - area = 0 - for i in range(n_corners): - j = (i + 1) % n_corners # next corner - area += X[i] * Y[j] - X[j] * Y[i] - area = abs(area) / 2 - return area + is_picked = is_loc_at(x, y, locs, r) + picked_locs = locs[is_picked] + return picked_locs @_numba.jit(nopython=True) @@ -302,13 +403,52 @@ def check_if_in_rectangle(x, y, X, Y): def locs_in_rectangle(locs, X, Y): + """Returns localizations in rectangle defined by corners (X, Y). + + Parameters + ---------- + locs : numpy.recarray + Localizations list. + X : list + x-coordinates of rectangle corners. + Y : list + y-coordinates of rectangle corners. + + Returns + ------- + picked_locs : numpy.recarray + Localizations in rectangle. + """ + is_in_rectangle = check_if_in_rectangle( locs.x, locs.y, _np.array(X), _np.array(Y) ) - return locs[is_in_rectangle] + picked_locs = locs[is_in_rectangle] + return picked_locs def minimize_shifts(shifts_x, shifts_y, shifts_z=None): + """Minimizes shifts in x, y, and z directions. Used for drift correction. + + Parameters + ---------- + shifts_x : numpy.2darray + Shifts in x direction. + shifts_y : numpy.2darray + Shifts in y direction. + shifts_z : numpy.2darray (default=None) + Shifts in z direction. + + Returns + ------- + shift_y : numpy.1darray + Minimized shifts in y direction. + shift_x : numpy.1darray + Minimized shifts in x direction. + shift_z : numpy.1darray (optional) + Minimized shifts in z direction if shifts_z is not None. + """ + n_channels = shifts_x.shape[0] n_pairs = int(n_channels * (n_channels - 1) / 2) n_dims = 2 if shifts_z is None else 3 @@ -334,11 +474,29 @@ def minimize_shifts(shifts_x, shifts_y, shifts_z=None): def n_futures_done(futures): + """Returns the number of finished futures, used in multiprocessing.""" + return sum([_.done() for _ in futures]) def remove_from_rec(rec_array, name): - return _drop_fields(rec_array, name, usemask=False, asrecarray=True) + """Removes a column from the existing np.recarray. + + Parameters + ---------- + rec_array : np.rec.array + Recarray from which the column is removed. + name : str + Name of the column to be removed. + + Returns + ------- + rec_array : np.rec.array + Recarray without the column. + """ + + rec_array = _drop_fields(rec_array, name, usemask=False, asrecarray=True) + return rec_array def locs_glob_map(func, pattern, args=[], kwargs={}, extension=""): @@ -359,3 +517,141 @@ def locs_glob_map(func, pattern, args=[], kwargs={}, extension=""): out_path = base + "_" + extension + ".hdf5" locs, info = result _io.save_locs(out_path, locs, info) + + +def get_pick_polygon_corners(pick): + """Returns X and Y coordinates of a pick polygon. + + Returns None, None if the pick is not a closed polygon.""" + + if len(pick) < 3 or pick[0] != pick[-1]: + return None, None + else: + X = [_[0] for _ in pick] + Y = [_[1] for _ in pick] + return X, Y + + +def get_pick_rectangle_corners(start_x, start_y, end_x, end_y, width): + """Finds the positions of corners of a rectangular pick. + Rectangular pick is defined by: + [(start_x, start_y), (end_x, end_y)] + and its width. (all values in camera pixels) + + Returns + ------- + corners : tuple + Contains corners' x and y coordinates in two lists + """ + + if end_x == start_x: + alpha = _np.pi / 2 + else: + alpha = _np.arctan((end_y - start_y) / (end_x - start_x)) + dx = width * _np.sin(alpha) / 2 + dy = width * _np.cos(alpha) / 2 + x1 = float(start_x - dx) + x2 = float(start_x + dx) + x4 = float(end_x - dx) + x3 = float(end_x + dx) + y1 = float(start_y + dy) + y2 = float(start_y - dy) + y4 = float(end_y + dy) + y3 = float(end_y - dy) + corners = ([x1, x2, x3, x4], [y1, y2, y3, y4]) + return corners + + +def pick_areas_circle(picks, r): + """Returns pick areas for each pick in picks. + + Parameters + ---------- + picks : list + List of picks, each pick is a list of x and y coordinates. + r : float + Pick radius. + + Returns + ------- + areas : np.1darray + Pick areas, same units as r. + """ + + areas = _np.ones(len(picks)) * _np.pi * r**2 + return areas + + +def polygon_area(X, Y): + """Finds the area of a polygon defined by corners X and Y. + + Parameters + ---------- + X : numpy.1darray + x-coordinates of the polygon corners. + Y : numpy.1darray + y-coordinates of the polygon corners. + + Returns + ------- + area : float + Area of the polygon. + """ + + n_corners = len(X) + area = 0 + for i in range(n_corners): + j = (i + 1) % n_corners # next corner + area += X[i] * Y[j] - X[j] * Y[i] + area = abs(area) / 2 + return area + + +def pick_areas_polygon(picks): + """Returns pick areas for each pick in picks. + + Parameters + ---------- + picks : list + List of picks, each pick is a list of coordinates of the + polygon corners. + + Returns + ------- + areas : np.1darray + Pick areas. + """ + + areas = _np.zeros(len(picks)) + for i, pick in enumerate(picks): + if len(pick) < 3 or pick[0] != pick[-1]: # not a closed polygon + areas[i] = 0 + continue + X, Y = get_pick_polygon_corners(pick) + areas[i] = polygon_area(X, Y) + areas = areas[areas > 0] # remove open polygons + return areas + + +def pick_areas_rectangle(picks, w): + """Returns pick areas for each pick in picks. + + Parameters + ---------- + picks : list + List of picks, each pick is a list of coordinates of the + rectangle corners. + w : float + Pick width. + + Returns + ------- + areas : np.1darray + Pick areas, same units as w. + """ + + areas = _np.zeros(len(picks)) + for i, pick in enumerate(picks): + (xs, ys), (xe, ye) = pick + areas[i] = w * _np.sqrt((xe - xs) ** 2 + (ye - ys) ** 2) + return areas \ No newline at end of file diff --git a/picasso/postprocess.py b/picasso/postprocess.py index ceee4c77..af74d7d9 100644 --- a/picasso/postprocess.py +++ b/picasso/postprocess.py @@ -94,6 +94,7 @@ def get_block_locs_at(x, y, index_blocks): indices = list(_itertools.chain(*indices)) return locs[indices] + @_numba.jit(nopython=True, nogil=True) def _fill_index_blocks( block_starts, block_ends, x_index, y_index @@ -107,6 +108,7 @@ def _fill_index_blocks( block_starts, block_ends, N, x_index, y_index, i, j, k ) + @_numba.jit(nopython=True, nogil=True) def _fill_index_block(block_starts, block_ends, N, x_index, y_index, i, j, k): block_starts[i, j] = k @@ -115,6 +117,144 @@ def _fill_index_block(block_starts, block_ends, N, x_index, y_index, i, j, k): block_ends[i, j] = k return k + +def picked_locs( + locs, info, picks, pick_shape, pick_size=None, add_group=True, callback=None +): + """Finds picked localizations. + + Parameters + ---------- + locs : np.recarray + Localization list. + info : list of dicts + Metadata of the localizations list. + picks : list + List of picks. + pick_shape : {'Circle', 'Rectangle', 'Polygon'} + Shape of the pick. + pick_size : float (default=None) + Size of the pick. Radius for the circles, width for the + rectangles, None for the polygons. + add_group : boolean (default=True) + True if group id should be added to locs. Each pick will be + assigned a different id. + callback : function (default=None) + Function to display progress. If "console", tqdm is used to + display the progress. If None, no progress is displayed. + + Returns + ------- + picked_locs : list of np.recarrays + List of np.recarrays, each containing locs from one pick. + """ + + if len(picks): + picked_locs = [] + if callback == "console": + progress = _tqdm(range(len(picks)), desc="Picking locs", unit="pick") + + if pick_shape == "Circle": + index_blocks = get_index_blocks(locs, info, pick_size) + for i, pick in enumerate(picks): + x, y = pick + block_locs = get_block_locs_at( + x, y, index_blocks + ) + group_locs = _lib.locs_at(x, y, block_locs, pick_size) + if add_group: + group = i * _np.ones(len(group_locs), dtype=_np.int32) + group_locs = _lib.append_to_rec( + group_locs, group, "group" + ) + group_locs.sort(kind="mergesort", order="frame") + picked_locs.append(group_locs) + + if callback == "console": + progress.update(1) + elif callback is not None: + callback(i + 1) + + elif pick_shape == "Rectangle": + for i, pick in enumerate(picks): + (xs, ys), (xe, ye) = pick + X, Y = _lib.get_pick_rectangle_corners( + xs, ys, xe, ye, pick_size + ) + x_min = min(X) + x_max = max(X) + y_min = min(Y) + y_max = max(Y) + group_locs = locs[locs.x > x_min] + group_locs = group_locs[group_locs.x < x_max] + group_locs = group_locs[group_locs.y > y_min] + group_locs = group_locs[group_locs.y < y_max] + group_locs = _lib.locs_in_rectangle(group_locs, X, Y) + # store rotated coordinates in x_rot and y_rot + angle = 0.5 *_np.pi - _np.arctan2((ye - ys), (xe - xs)) + x_shifted = group_locs.x - xs + y_shifted = group_locs.y - ys + x_pick_rot = x_shifted * _np.cos( + angle + ) - y_shifted * _np.sin(angle) + y_pick_rot = x_shifted * _np.sin( + angle + ) + y_shifted * _np.cos(angle) + group_locs = _lib.append_to_rec( + group_locs, x_pick_rot, "x_pick_rot" + ) + group_locs = _lib.append_to_rec( + group_locs, y_pick_rot, "y_pick_rot" + ) + if add_group: + group = i * _np.ones(len(group_locs), dtype=_np.int32) + group_locs = _lib.append_to_rec( + group_locs, group, "group" + ) + group_locs.sort(kind="mergesort", order="frame") + picked_locs.append(group_locs) + + if callback == "console": + progress.update(1) + elif callback is not None: + callback(i + 1) + + elif pick_shape == "Polygon": + for i, pick in enumerate(picks): + X, Y = _lib.get_pick_polygon_corners(pick) + if X is None: + if callback == "console": + progress.update(1) + elif callback is not None: + callback(i + 1) + continue + group_locs = locs[locs.x > min(X)] + group_locs = group_locs[group_locs.x < max(X)] + group_locs = group_locs[group_locs.y > min(Y)] + group_locs = group_locs[group_locs.y < max(Y)] + group_locs = _lib.locs_in_polygon(group_locs, X, Y) + if add_group: + group = i * _np.ones(len(group_locs), dtype=_np.int32) + group_locs = _lib.append_to_rec( + group_locs, group, "group" + ) + group_locs.sort(kind="mergesort", order="frame") + picked_locs.append(group_locs) + + if callback == "console": + progress.update(1) + elif callback is not None: + callback(i + 1) + + else: + raise ValueError( + "Invalid pick shape. Please choose from 'Circle', 'Rectangle', " + "'Polygon'." + ) + + return picked_locs + + @_numba.jit(nopython=True, nogil=True, cache=True) def pick_similar( x, y_shift, y_base, @@ -189,6 +329,7 @@ def pick_similar( ) return x_similar, y_similar + @_numba.jit(nopython=True, nogil=True) def _n_block_locs_at(x_range, y_range, K, L, block_starts, block_ends, cache=True): step = 0 @@ -203,6 +344,7 @@ def _n_block_locs_at(x_range, y_range, K, L, block_starts, block_ends, cache=Tru n_block_locs += _np.uint32(block_ends[k][l] - block_starts[k][l]) return n_block_locs + @_numba.jit(nopython=True, nogil=True, cache=True) def _get_block_locs_at( x_range, y_range, locs_xy, @@ -226,6 +368,7 @@ def _get_block_locs_at( )) return locs_xy[:, indices] + @_numba.jit(nopython=True, nogil=True, cache=True) def _locs_at(x, y, locs_xy, r): dx = locs_xy[0] - x @@ -234,12 +377,14 @@ def _locs_at(x, y, locs_xy, r): is_picked = dx ** 2 + dy ** 2 < r2 return locs_xy[:, is_picked] + @_numba.jit(nopython=True, nogil=True) def _rmsd_at_com(locs_xy): com_x = _np.mean(locs_xy[0]) com_y = _np.mean(locs_xy[1]) return _np.sqrt(_np.mean((locs_xy[0] - com_x) ** 2 + (locs_xy[1] - com_y) ** 2)) + @_numba.jit(nopython=True, nogil=True) def _distance_histogram( locs, diff --git a/picasso/zfit.py b/picasso/zfit.py index e61f3bd3..082a4d2b 100644 --- a/picasso/zfit.py +++ b/picasso/zfit.py @@ -61,6 +61,7 @@ def calibrate_z(locs, info, d, magnification_factor, path=None): "Y Coefficients": [float(_) for _ in cy], "Number of frames": int(n_frames), "Step size in nm": float(d), + "Magnification factor": float(magnification_factor), } if path is not None: with open(path, "w") as f: