From c1dff68288ecfefd5a3c9371b0f03f5381f867cf Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Thu, 6 Jun 2024 17:33:40 +0200 Subject: [PATCH 01/16] modify pypi testing action FIX --- .github/workflows/pypi_release.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/pypi_release.yml b/.github/workflows/pypi_release.yml index 0596e4eb..f9c9c34f 100644 --- a/.github/workflows/pypi_release.yml +++ b/.github/workflows/pypi_release.yml @@ -50,6 +50,7 @@ jobs: - uses: conda-incubator/setup-miniconda@v2 with: auto-update-conda: true + miniconda-version: "latest" python-version: ${{ matrix.python-version }} - name: Conda info shell: bash -l {0} From 70d6d9accda65187fd967413214c58231e690f97 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 28 Jun 2024 15:10:56 +0200 Subject: [PATCH 02/16] fix clusterfile (CMD) --- changelog.rst | 6 +++++- picasso/io.py | 5 ++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/changelog.rst b/changelog.rst index 16e9f796..e9788a85 100644 --- a/changelog.rst +++ b/changelog.rst @@ -1,7 +1,11 @@ Changelog ========= -Last change: 06-JUN-2024 MTS +Last change: 28-JUN-2024 MTS + +0.7.0 +----- +- CMD ``clusterfile`` fixed 0.6.9 - 0.6.11 -------------- diff --git a/picasso/io.py b/picasso/io.py index b193a025..3c052b34 100644 --- a/picasso/io.py +++ b/picasso/io.py @@ -1249,7 +1249,10 @@ def load_locs(path, qt_parent=None): def load_clusters(path, qt_parent=None): with _h5py.File(path, "r") as cluster_file: - clusters = cluster_file["clusters"][...] + try: + clusters = cluster_file["clusters"][...] + except KeyError: + clusters = cluster_file["locs"][...] clusters = _np.rec.array( clusters, dtype=clusters.dtype ) # Convert to rec array with fields as attributes From d52fd18e07bfa1e5d623b8a98986661df38be0f3 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 5 Jul 2024 13:51:42 +0200 Subject: [PATCH 03/16] implementation of AIM in Picasso --- changelog.rst | 1 + picasso/aim.py | 705 ++++++++++++++++++++++++++++++++++++++++ picasso/gui/render.py | 475 ++++++++++++++++----------- picasso/imageprocess.py | 2 +- picasso/lib.py | 7 + readme.rst | 8 +- 6 files changed, 998 insertions(+), 200 deletions(-) create mode 100644 picasso/aim.py diff --git a/changelog.rst b/changelog.rst index e9788a85..7a2a3d7d 100644 --- a/changelog.rst +++ b/changelog.rst @@ -5,6 +5,7 @@ Last change: 28-JUN-2024 MTS 0.7.0 ----- +- Adaptive Intersection Maximization (AIM, doi: 10.1038/s41592-022-01307-0) implemented - CMD ``clusterfile`` fixed 0.6.9 - 0.6.11 diff --git a/picasso/aim.py b/picasso/aim.py new file mode 100644 index 00000000..650969ef --- /dev/null +++ b/picasso/aim.py @@ -0,0 +1,705 @@ +""" + picasso.aim + ~~~~~~~~~~~ + + Picasso implementation of Adaptive Intersection Maximization (AIM) + for fast undrifting in 2D and 3D. + + Adapted from: Ma, H., et al. Science Advances. 2024. + + :author: Hongqiang Ma, Maomao Chen, Phuong Nguyen, Yang Liu, + Rafal Kowalewski, 2024 + :copyright: Copyright (c) 2016-2024 Jungmann Lab, MPI of Biochemistry +""" + +from concurrent.futures import ThreadPoolExecutor as _ThreadPoolExecutor + +import numpy as _np +from scipy.interpolate import InterpolatedUnivariateSpline as \ + _InterpolatedUnivariateSpline +from tqdm import tqdm as _tqdm + +#TODO: add progress dialog, also tqdm + + +def intersect1d(a, b): + """Slightly faster implementation of _np.intersect1d without + unnecessary checks, etc. + + Finds the indices of common elements in two 1D arrays (a and b). + Both a and b are assumed to be sorted and contain only unique + values. + + Parameters + ---------- + a : _np.array + 1D array of integers. + b : _np.array + 1D array of integers. + + Returns + ------- + a_indices : _np.array + Indices of common elements in a. + b_indices : _np.array + Indices of common elements in b. + """ + + aux = _np.concatenate((a, b)) + aux_sort_indices = _np.argsort(aux, kind='mergesort') + aux = aux[aux_sort_indices] + + mask = aux[1:] == aux[:-1] + a_indices = aux_sort_indices[:-1][mask] + b_indices = aux_sort_indices[1:][mask] - a.size + + return a_indices, b_indices + + +def count_intersections(l0_coords, l0_counts, l1_coords, l1_counts): + """Counts the number of intersected localizations between the two + datasets. We assume that the intersection distance is 1 and since + the coordinates are expressed in the units of intersection distance, + we require the coordinates to be exactly the same to count as + intersection. Also, coordinates are converted to 1D arrays + (x + y * width). + + Parameters + ---------- + l0_coords : _np.array + Unique coordinates of the reference localizations. + l0_counts : _np.array + Counts of the unique values of reference localizations. + l1_coords : _np.array + Unique coordinates of the target localizations. + l1_counts : _np.array + Counts of the unique values of target localizations. + + Returns + ------- + n_intersections : int + Number of intersections. + """ + + # indices of common elements + idx0, idx1 = intersect1d(l0_coords, l1_coords) + # extract the counts of these elements + l0_counts_subset = l0_counts[idx0] + l1_counts_subset = l1_counts[idx1] + # for each overlapping coordinate, take the minimum count from l0 + # and l1, sum up across all overlapping coordinates + n_intersections = _np.sum( + _np.minimum(l0_counts_subset, l1_counts_subset) + ) + return n_intersections + + +def run_intersections( + l0_coords, l0_counts, l1_coords, l1_counts, shifts_xy, box +): + """Run intersection counting across the local search region. Returns + the 2D array with number of intersections across the local search + region. + + Parameters + ---------- + l0_coords : _np.array + Unique coordinates of the reference localizations. + l0_counts : _np.array + Counts of the reference localizations. + l1_coords : _np.array + Unique coordinates of the target localizations. + l1_counts : _np.array + Counts of the target localizations. + shifts_xy : _np.array + 1D array with x and y shifts. + box : int + Side length of the local search region. + + Returns + ------- + roi_cc : _np.2darray + 2D array with number of intersections across the local search + region. + """ + + # create the 2D array with shifts + roi_cc = _np.zeros(shifts_xy.shape, dtype=_np.int32) + # shift target coordinates + l1_coords_shifted = l1_coords[:, _np.newaxis] + shifts_xy + # go through each element in the local search region + for i in range(len(shifts_xy)): + n_intersections = count_intersections( + l0_coords, l0_counts, l1_coords_shifted[:, i], l1_counts + ) + roi_cc[i] = n_intersections + return roi_cc.reshape(box, box) + + +def run_intersections_multithread( + l0_coords, l0_counts, l1_coords, l1_counts, shifts_xy, box +): + """Run intersection counting across the local search region. Returns + the 2D array with number of intersections across the local search + region. Uses multithreading. + + Parameters + ---------- + l0_coords : _np.array + Unique coordinates of the reference localizations. + l0_counts : _np.array + Counts of the reference localizations. + l1_coords : _np.array + Unique coordinates of the target localizations. + l1_counts : _np.array + Counts of the target localizations. + shifts_xy : _np.array + 1D array with x and y shifts. + box : int + Side length of the local search region. + + Returns + ------- + roi_cc : _np.2darray + 2D array with number of intersections across the local search + region. + """ + + # shift target coordinates + l1_coords_shifted = l1_coords[:, _np.newaxis] + shifts_xy + # run multiple threads + n_workers = len(shifts_xy) + executor = _ThreadPoolExecutor(n_workers) + f = [ + executor.submit( + count_intersections, + l0_coords, + l0_counts, + l1_coords_shifted[:, i], + l1_counts, + ) + for i in range(len(shifts_xy)) + ] + executor.shutdown(wait=True) + if box == 1: # z intersection only, for z undrifting + roi_cc = _np.array([_.result() for _ in f]) + else: # 2D intersection + roi_cc = _np.array([_.result() for _ in f]).reshape(box, box) + return roi_cc + + +def point_intersect_2d( + l0_coords, l0_counts, x1, y1, intersect_d, width_units, shifts_xy, box +): + """Converts target coordinates into a 1D array in units of + intersect_d and counts the number of intersections in the local + search region. + + Parameters + ---------- + l0_coords : _np.array + Unique values of the reference localizations. + l0_counts : _np.array + Counts of the unique values of reference localizations. + x1 : _np.array + x-coordinates of the target (currently undrifted) localizations. + y1 : _np.array + y-coordinates of the target (currently undrifted) localizations. + intersect_d : float + Intersect distance in camera pixels. + width_units : int + Width of the camera image in units of intersect_d. + shifts_xy : _np.array + 1D array with x and y shifts. + box : int + Final side length of the local search region. + + Returns + ------- + roi_cc : _np.2darray + 2D array with numbers of intersections in the local search + region. + """ + + # convert target coordinates to a 1D array in intersect_d units + x1_units = _np.round(x1 / intersect_d) + y1_units = _np.round(y1 / intersect_d) + l1 = _np.int32(x1_units + y1_units * width_units) # 1d list + # get unique values and counts of the target localizations + l1_coords, l1_counts = _np.unique(l1, return_counts=True) + # run the intersections counting + roi_cc = run_intersections_multithread( + l0_coords, l0_counts, l1_coords, l1_counts, shifts_xy, box + ) + return roi_cc + + +def point_intersect_3d( + l0_coords, l0_counts, x1, y1, z1, + intersect_d, width_units, height_units, shifts_z, +): + """Converts target coordinates into a 1D array in units of + intersect_d and counts the number of intersections in the local + search region. + + Parameters + ---------- + l0_coords : _np.array + Unique values of the reference localizations. + l0_counts : _np.array + Counts of the unique values of reference localizations. + x1 : _np.array + x-coordinates of the target (currently undrifted) localizations. + y1 : _np.array + y-coordinates of the target (currently undrifted) localizations. + z1 : _np.array + z-coordinates of the target (currently undrifted) localizations. + intersect_d : float + Intersect distance in camera pixels. + width_units : int + Width of the camera image in units of intersect_d. + height_units : int + Height of the camera image in units of intersect_d. + shifts_z : _np.array + 1D array with z shifts. + + Returns + ------- + roi_cc : _np.2darray + 2D array with numbers of intersections in the local search + region. + """ + + # convert target coordinates to a 1D array in intersect_d units + x1_units = _np.round(x1 / intersect_d) + y1_units = _np.round(y1 / intersect_d) + z1_units = _np.round(z1 / intersect_d) + l1 = _np.int32( + x1_units + + y1_units * width_units + + z1_units * width_units * height_units + ) # 1d list + # get unique values and counts of the target localizations + l1_coords, l1_counts = _np.unique(l1, return_counts=True) + # run the intersections counting + roi_cc = run_intersections_multithread( + l0_coords, l0_counts, l1_coords, l1_counts, shifts_z, 1 + ) + return roi_cc + + +def get_fft_peak(roi_cc, roi_size): + """Estimate the precise sub-pixel position of the peak of roi_cc + with FFT. + + Parameters + ---------- + roi_cc : _np.2darray + 2D array with numbers of intersections in the local search region. + roi_size : int + Size of the local search region. + + Returns + ------- + px : float + Estimated x-coordinate of the peak. + py : float + Estimated y-coordinate of the peak. + """ + + fft_values = _np.fft.fft2(roi_cc.T) + ang_x = _np.angle(fft_values[0, 1]) + ang_x = ang_x - 2 * _np.pi * (ang_x > 0) # normalize + px = ( + _np.abs(ang_x) / (2 * _np.pi / roi_cc.shape[0]) + - (roi_cc.shape[0] - 1) / 2 + ) # peak in x + px *= roi_size / roi_cc.shape[0] # convert to intersect_d units + ang_y = _np.angle(fft_values[1, 0]) + ang_y = ang_y - 2 * _np.pi * (ang_y > 0) # normalize + py = ( + _np.abs(ang_y) / (2 * _np.pi / roi_cc.shape[1]) + - (roi_cc.shape[1] - 1) / 2 + ) # peak in y + py *= roi_size / roi_cc.shape[1] # convert to intersect_d units + return px, py + + +def get_fft_peak_z(roi_cc, roi_size): + """Estimate the precise sub-pixel position of the peak of 1D roi_cc. + + Parameters + ---------- + roi_cc : _np.array + 1D array with numbers of intersections in the local search + region. + roi_size : int + Size of the local search region. + + Returns + ------- + pz : float + Estimated z-coordinate of the peak. + """ + + fft_values = _np.fft.fft(roi_cc) + ang_z = _np.angle(fft_values[1]) + ang_z = ang_z - 2 * _np.pi * (ang_z > 0) # normalize + pz = ( + _np.abs(ang_z) / (2 * _np.pi / roi_cc.size) + - (roi_cc.size - 1) / 2 + ) # peak in z + pz *= roi_size / roi_cc.size # convert to intersect_d units + return pz + + +def intersection_max( + x, y, ref_x, ref_y, + frame, segmentation, intersect_d, roi_r, width, progress=None, +): + """Maximize intersection (undrift) for 2D localizations. + + Parameters + ---------- + x : _np.array + x-coordinates of the localizations. + y : _np.array + y-coordinates of the localizations. + ref_x_list : _np.array + x-coordinates of the reference localizations. + ref_y_list : _np.array + y-coordinates of the reference localizations. + frame : _np.array + Frame indices of localizations. + segmentation : int + Time interval for drift tracking, unit: frames. + intersect_d : float + Intersect distance in camera pixels. + roi_r : float + Radius of the local search region in camera pixels. Should be + higher than the maximum expected drift within one segment. + width : int + Width of the camera image in camera pixels. + progress : picasso.lib.ProgressDialog (default=None) + Progress dialog. If None, progress is displayed with tqdm. + + Returns + ------- + x_pdc : _np.array + Undrifted x-coordinates. + y_pdc : _np.array + Undrifted y-coordinates. + drift_x : _np.array + Drift in x-direction. + drift_y : _np.array + Drift in y-direction. + """ + + # number of segments + n_segments = int(_np.ceil(frame.max() / segmentation)) + rel_drift_x = 0 # adaptive drift (updated at each interval) + rel_drift_y = 0 + + # drift in x and y + drift_x = _np.zeros(n_segments) + drift_y = _np.zeros(n_segments) + + # find shifts for the local search region (in units of intersect_d) + roi_units = int(_np.ceil(roi_r / intersect_d)) + steps = _np.arange(-roi_units, roi_units + 1, 1) + box = len(steps) + shifts_xy = _np.zeros((box, box), dtype=_np.int32) + width_units = width / intersect_d + for i, shift_x in enumerate(steps): + for j, shift_y in enumerate(steps): + shifts_xy[i, j] = shift_x + shift_y * width_units + shifts_xy = shifts_xy.reshape(box**2) + + # convert reference to a 1D array in units of intersect_d and find + # unique values and counts + x0_units = _np.round(ref_x / intersect_d) + y0_units = _np.round(ref_y / intersect_d) + l0 = _np.int32(x0_units + y0_units * width_units) # 1d list + l0_coords, l0_counts = _np.unique(l0, return_counts=True) + + # initialize progress such that if GUI is used, tqdm is omitted + if progress is not None: + iterator = range(1, n_segments) + else: + iterator = _tqdm(n_segments, desc="Undrifting z", unit="segment") + + # run across each segment + for s in iterator: + # get the target localizations within the current segment + min_frame_idx = frame > s * segmentation + max_frame_idx = frame <= (s + 1) * segmentation + x1 = x[min_frame_idx & max_frame_idx] + y1 = y[min_frame_idx & max_frame_idx] + + # undrifting from the previous round + x1 += rel_drift_x + y1 += rel_drift_y + + # count the number of intersected localizations + roi_cc = point_intersect_2d( + l0_coords, l0_counts, x1, y1, + intersect_d, width_units, shifts_xy, box, + ) + + # estimate the precise sub-pixel position of the peak of roi_cc + # with FFT + px, py = get_fft_peak(roi_cc, 2 * roi_r) + + # update the relative drift reference for the subsequent + # segmented subset (interval) and save the drifts + rel_drift_x += px + rel_drift_y += py + drift_x[s] = -rel_drift_x + drift_y[s] = -rel_drift_y + + # update progress + if progress is not None: + progress.set_value(s) + else: + iterator.update(s - iterator.n) + + # interpolate the drifts (cubic spline) for all frames + n_frames = n_segments * segmentation + drift_track_points = _np.linspace(0, n_frames, n_segments+1) + t = (drift_track_points[1:] + drift_track_points[:-1]) / 2 + drift_x_pol = _InterpolatedUnivariateSpline(t, drift_x, k=3) + drift_y_pol = _InterpolatedUnivariateSpline(t, drift_y, k=3) + t_inter = _np.arange(n_frames) + 1 + drift_x = drift_x_pol(t_inter) + drift_y = drift_y_pol(t_inter) + + # undrift the localizations + x_pdc = x - drift_x[frame-1] + y_pdc = y - drift_y[frame-1] + + return x_pdc, y_pdc, drift_x, drift_y + + +def intersection_max_z( + x, y, z, ref_x, ref_y, ref_z, + frame, segmentation, intersect_d, roi_r, width, height, pixelsize, + progress=None, +): + """Maximize intersection (undrift) for 3D localizations. Assumes + that x and y coordinates were already undrifted. x and y are in + units of camera pixels, z is in nm. + + See intersection_max for more details.""" + + # convert z to camera pixels + z = z.copy() / pixelsize + ref_z = ref_z.copy() / pixelsize #TODO: remember to convert back to nm, also for drift_z + + # number of segments + n_segments = int(_np.ceil(frame.max() / segmentation)) + rel_drift_z = 0 # adaptive drift (updated at each interval) + + # drift in z + drift_z = _np.zeros(n_segments) + + # find shifts for the local search region (in units of intersect_d) + roi_units = int(_np.ceil(roi_r / intersect_d)) + steps = _np.arange(-roi_units, roi_units + 1, 1) + box = len(steps) + width_units = width / intersect_d + height_units = height / intersect_d + shifts_z = steps.astype(_np.int32) * width_units * height_units + + # convert reference to a 1D array in units of intersect_d and find + # unique values and counts + x0_units = _np.round(ref_x / intersect_d) + y0_units = _np.round(ref_y / intersect_d) + z0_units = _np.round(ref_z / intersect_d) + l0 = _np.int32( + x0_units + + y0_units * width_units + + z0_units * width_units * height_units + ) # 1d list + l0_coords, l0_counts = _np.unique(l0, return_counts=True) + + # initialize progress such that if GUI is used, tqdm is omitted + if progress is not None: + iterator = range(1, n_segments) + else: + iterator = _tqdm(n_segments, desc="Undrifting z", unit="segment") + + # run across each segment + for s in iterator: + # get the target localizations within the current segment + min_frame_idx = frame > s * segmentation + max_frame_idx = frame <= (s + 1) * segmentation + x1 = x[min_frame_idx & max_frame_idx] + y1 = y[min_frame_idx & max_frame_idx] + z1 = z[min_frame_idx & max_frame_idx] + + # undrifting from the previous round + z1 += rel_drift_z + + # count the number of intersected localizations + roi_cc = point_intersect_3d( + l0_coords, l0_counts, x1, y1, z1, + intersect_d, width_units, height_units, shifts_z, + ) + + # estimate the precise sub-pixel position of the peak of roi_cc + # with FFT + pz = get_fft_peak_z(roi_cc, 2 * roi_r) + + # update the relative drift reference for the subsequent + # segmented subset (interval) and save the drifts + rel_drift_z += pz + drift_z[s] = -rel_drift_z + + # update progress + if progress is not None: + progress.set_value(s) + else: + iterator.update(s - iterator.n) + + + # interpolate the drifts (cubic spline) for all frames + n_frames = n_segments * segmentation + drift_track_points = _np.linspace(0, n_frames, n_segments+1) + t = (drift_track_points[1:] + drift_track_points[:-1]) / 2 + drift_z_pol = _InterpolatedUnivariateSpline(t, drift_z, k=3) + t_inter = _np.arange(n_frames) + 1 + drift_z = drift_z_pol(t_inter) + + # undrift the localizations + z_pdc = z - drift_z[frame-1] + + # convert back to nm + z_pdc *= pixelsize + drift_z *= pixelsize + + return z_pdc, drift_z + + +def aim( + locs, + info, + segmentation=100, + intersect_d=20/130, + roi_r=60/130, + progress=None, +): + """Apply AIM algorithm to the localizations. + + Parameters + ---------- + locs : _np.rec.array + Localizations list to be undrifted. + info : list of dicts + Localizations list's metadata. + intersect_d : float + Intersect distance in camera pixels. + segmentation : int + Time interval for drift tracking, unit: frames. + roi_r : float + Radius of the local search region in camera pixels. Should be + larger than the maximum expected drift within segmentation. + progress : picasso.lib.ProgressDialog (default=None) + Progress dialog. If None, progress is displayed with tqdm. + + Returns + ------- + locs : _np.rec.array + Undrifted localizations. + drift : _np.rec.array + Drift in x and y directions (and z if applicable). + """ + + # extract metadata + width = info[0]["Width"] + height = info[0]["Height"] + pixelsize = info[1]["Pixelsize"] + + # frames should start at 1 + frame = locs["frame"] + 1 + n_frames = _np.max(frame) + + # group frames into track intervals and correct for rounding errors + n_segments = _np.floor(n_frames / segmentation).astype(int) + n_frames = n_segments * segmentation + frame[frame > n_frames] = n_frames # cap frame number + + # get the reference localizations (first interval) + ref_x = locs["x"][frame <= segmentation] + ref_y = locs["y"][frame <= segmentation] + + ### RUN AIM TWICE ### + # the first run is with the first interval as reference + x_pdc, y_pdc, drift_x1, drift_y1 = intersection_max( + locs.x, locs.y, ref_x, ref_y, + frame, segmentation, intersect_d, roi_r, width, + progress=progress, + ) + # the second run is with the entire dataset as reference + if progress is not None: + progress.zero_progress(description="Undrifting by AIM (2/2)") + x_pdc, y_pdc, drift_x2, drift_y2 = intersection_max( + x_pdc, y_pdc, x_pdc, y_pdc, + frame, segmentation, intersect_d, roi_r, width, + progress=progress, + ) + + # add the drifts together from the two rounds + drift_x = drift_x1 + drift_x2 + drift_y = drift_y1 + drift_y2 + + # shift the drifts by the mean value (like in Picasso) + drift_x -= _np.mean(drift_x) + drift_y -= _np.mean(drift_y) + + # combine to Picasso format + drift = _np.rec.array((drift_x, drift_y), dtype=[("x", "f"), ("y", "f")]) + + # 3D undrifting + if hasattr(locs, "z"): + if progress is not None: + progress.zero_progress(description="Undrifting z (1/2)") + ref_x = x_pdc[frame <= segmentation] + ref_y = y_pdc[frame <= segmentation] + ref_z = locs.z[frame <= segmentation] + z_pdc, drift_z1 = intersection_max_z( + x_pdc, y_pdc, locs.z, ref_x, ref_y, ref_z, + frame, segmentation, intersect_d, roi_r, width, height, pixelsize, + progress=progress, + ) + if progress is not None: + progress.zero_progress(description="Undrifting z (2/2)") + z_pdc, drift_z2 = intersection_max_z( + x_pdc, y_pdc, z_pdc, x_pdc, y_pdc, z_pdc, + frame, segmentation, intersect_d, roi_r, width, height, pixelsize, + progress=progress, + ) + drift_z = drift_z1 + drift_z2 + drift_z -= _np.mean(drift_z) + drift = _np.rec.array( + (drift_x, drift_y, drift_z), + dtype=[("x", "f"), ("y", "f"), ("z", "f")] + ) + + # apply the drift to localizations + locs["x"] = x_pdc + locs["y"] = y_pdc + if hasattr(locs, "z"): + locs["z"] = z_pdc + + new_info = { + "Undrifted by": "AIM", + "Intersect distance (cam. pixels)": intersect_d, + "Segmentation": segmentation, + "Search regions radius (cam. pixels)": roi_r, + } + new_info = info + [new_info] + + if progress is not None: + progress.close() + + return locs, new_info, drift \ No newline at end of file diff --git a/picasso/gui/render.py b/picasso/gui/render.py index d2419a11..4600c954 100644 --- a/picasso/gui/render.py +++ b/picasso/gui/render.py @@ -2,8 +2,8 @@ gui/render ~~~~~~~~~~~~~~~~~~~~ Graphical user interface for rendering localization images - :author: Joerg Schnitzbauer & Maximilian Strauss - & Rafal Kowalewski, 2017-2022 + :author: Joerg Schnitzbauer, Maximilian Strauss, Rafal Kowalewski, + 2017-2022 :copyright: Copyright (c) 2017 Jungmann Lab, MPI of Biochemistry """ import os @@ -37,7 +37,8 @@ import colorsys -from .. import imageprocess, io, lib, postprocess, render, clusterer, __version__ +from .. import imageprocess, io, lib, postprocess, render, clusterer, aim, \ + __version__ from .rotation import RotationWindow # PyImarisWrite works on windows only @@ -1678,26 +1679,41 @@ def getParams(all_picked_locs, current, length, n_clusters, color_sys): clustered_locs, ) + @staticmethod + def getParams(parent=None): + """ + Creates the dialog and returns the requested values for + linking. + """ -class LinkDialog(QtWidgets.QDialog): - """ - A class to obtain inputs for linking localizations. + dialog = LinkDialog(parent) + result = dialog.exec_() + return ( + dialog.max_distance.value(), + dialog.max_dark_time.value(), + result == QtWidgets.QDialog.Accepted, + ) + + +class AIMDialog(QtWidgets.QDialog): + """Dialog to choose parameters for AIM undrifting. ... Attributes ---------- - max_dark_time : QDoubleSpinBox - contains the maximum gap between localizations (frames) to be - considered as belonging to the same group of linked locs - max_distance : QDoubleSpinBox - contains the maximum distance (pixels) between locs to be - considered as belonging to the same group of linked locs + intersect_d : QDoubleSpinBox + Contains the intersection distance in nm. + max_drift : QDoubleSpinBox + Contains the maximum drift within segmentation in nm. + segmentation : QSpinBox + Contains the lenght of temporal segments in units of frames. Methods ------- getParams(parent=None) - Creates the dialog and returns the requested values for linking + Creates the dialog and converts and returns the requested values + for AIM. """ def __init__(self, window): @@ -1706,19 +1722,33 @@ def __init__(self, window): self.setWindowTitle("Enter parameters") vbox = QtWidgets.QVBoxLayout(self) grid = QtWidgets.QGridLayout() - grid.addWidget(QtWidgets.QLabel("Max. distance (pixels):"), 0, 0) - self.max_distance = QtWidgets.QDoubleSpinBox() - self.max_distance.setRange(0, 1e6) - self.max_distance.setValue(1) - grid.addWidget(self.max_distance, 0, 1) - grid.addWidget(QtWidgets.QLabel("Max. transient dark frames:"), 1, 0) - self.max_dark_time = QtWidgets.QDoubleSpinBox() - self.max_dark_time.setRange(0, 1e9) - self.max_dark_time.setValue(1) - grid.addWidget(self.max_dark_time, 1, 1) + grid.addWidget(QtWidgets.QLabel("Segmentation:"), 0, 0) + self.segmentation = QtWidgets.QSpinBox() + self.segmentation.setRange(1, int(1e5)) + self.segmentation.setValue(100) + grid.addWidget(self.segmentation, 0, 1) + grid.addWidget(QtWidgets.QLabel("Intersection distance (nm):"), 1, 0) + self.intersect_d = QtWidgets.QDoubleSpinBox() + self.intersect_d.setRange(0.1, 1e6) + try: + default = 6 * float(window.info_dialog.fit_precision.text()) + except: + default = 20.0 + self.intersect_d.setValue(default) + self.intersect_d.setDecimals(1) + self.intersect_d.setSingleStep(1) + grid.addWidget(self.intersect_d, 1, 1) + grid.addWidget( + QtWidgets.QLabel("Max. drift in segmentation (nm):"), 2, 0 + ) + self.max_drift = QtWidgets.QDoubleSpinBox() + self.max_drift.setRange(0.1, 1e6) + self.max_drift.setValue(60.0) + self.max_drift.setDecimals(1) + self.max_drift.setSingleStep(1) + grid.addWidget(self.max_drift, 2, 1) vbox.addLayout(grid) - hbox = QtWidgets.QHBoxLayout() - vbox.addLayout(hbox) + # OK and Cancel buttons self.buttons = QtWidgets.QDialogButtonBox( QtWidgets.QDialogButtonBox.Ok | QtWidgets.QDialogButtonBox.Cancel, @@ -1731,18 +1761,18 @@ def __init__(self, window): @staticmethod def getParams(parent=None): - """ - Creates the dialog and returns the requested values for - linking. - """ + """Creates the dialog and converts and returns the requested + values for AIM.""" - dialog = LinkDialog(parent) + dialog = AIMDialog(parent) result = dialog.exec_() - return ( - dialog.max_distance.value(), - dialog.max_dark_time.value(), - result == QtWidgets.QDialog.Accepted, - ) + # convert intersect_d and max_drift to pixels + params = { + "segmentation": dialog.segmentation.value(), + "intersect_d": dialog.intersect_d.value(), + "roi_r": dialog.max_drift.value(), + } + return params, result == QtWidgets.QDialog.Accepted class DbscanDialog(QtWidgets.QDialog): @@ -1900,20 +1930,69 @@ def getParams(parent=None): dialog.save_centers.isChecked(), result == QtWidgets.QDialog.Accepted, ) + -class SMLMDialog3D(QtWidgets.QDialog): +class LinkDialog(QtWidgets.QDialog): """ - A class to obtain inputs for SMLM clusterer (3D). + A class to obtain inputs for linking localizations. ... Attributes ---------- - radius_xy : QDoubleSpinBox + max_dark_time : QDoubleSpinBox + contains the maximum gap between localizations (frames) to be + considered as belonging to the same group of linked locs + max_distance : QDoubleSpinBox + contains the maximum distance (pixels) between locs to be + considered as belonging to the same group of linked locs + + Methods + ------- + getParams(parent=None) + Creates the dialog and returns the requested values for linking + """ + + def __init__(self, window): + super().__init__(window) + self.window = window + self.setWindowTitle("Enter parameters") + vbox = QtWidgets.QVBoxLayout(self) + grid = QtWidgets.QGridLayout() + grid.addWidget(QtWidgets.QLabel("Max. distance (pixels):"), 0, 0) + self.max_distance = QtWidgets.QDoubleSpinBox() + self.max_distance.setRange(0, 1e6) + self.max_distance.setValue(1) + grid.addWidget(self.max_distance, 0, 1) + grid.addWidget(QtWidgets.QLabel("Max. transient dark frames:"), 1, 0) + self.max_dark_time = QtWidgets.QDoubleSpinBox() + self.max_dark_time.setRange(0, 1e9) + self.max_dark_time.setValue(1) + grid.addWidget(self.max_dark_time, 1, 1) + vbox.addLayout(grid) + hbox = QtWidgets.QHBoxLayout() + vbox.addLayout(hbox) + # OK and Cancel buttons + self.buttons = QtWidgets.QDialogButtonBox( + QtWidgets.QDialogButtonBox.Ok | QtWidgets.QDialogButtonBox.Cancel, + QtCore.Qt.Horizontal, + self, + ) + vbox.addWidget(self.buttons) + self.buttons.accepted.connect(self.accept) + self.buttons.rejected.connect(self.reject) + + +class SMLMDialog2D(QtWidgets.QDialog): + """ + A class to obtain inputs for SMLM clusterer (2D). + + ... + + Attributes + ---------- + radius : QDoubleSpinBox contains clustering radius in x and y directions - radius_z : QDoubleSpinBox - contains clustering radius in z direction - (typically =2*radius_xy) min_locs : QSpinBox contains minimum number of locs in cluster @@ -1923,43 +2002,36 @@ class SMLMDialog3D(QtWidgets.QDialog): Creates the dialog and returns the requested values for clustering """ - + def __init__(self, window): super().__init__(window) self.window = window - self.setWindowTitle("Enter parameters (3D)") + self.setWindowTitle("Enter parameters (2D)") vbox = QtWidgets.QVBoxLayout(self) grid = QtWidgets.QGridLayout() - # radius xy - grid.addWidget(QtWidgets.QLabel("Cluster radius xy (pixels):"), 0, 0) - self.radius_xy = QtWidgets.QDoubleSpinBox() - self.radius_xy.setRange(0.0001, 1e3) - self.radius_xy.setDecimals(4) - self.radius_xy.setValue(0.1) - grid.addWidget(self.radius_xy, 0, 1) - # radius z - grid.addWidget(QtWidgets.QLabel("Cluster radius z (pixels):"), 1, 0) - self.radius_z = QtWidgets.QDoubleSpinBox() - self.radius_z.setRange(0, 1e3) - self.radius_z.setDecimals(4) - self.radius_z.setValue(0.25) - grid.addWidget(self.radius_z, 1, 1) + # clustering radius + grid.addWidget(QtWidgets.QLabel("Cluster radius (pixels):"), 0, 0) + self.radius = QtWidgets.QDoubleSpinBox() + self.radius.setRange(0.0001, 1e3) + self.radius.setDecimals(4) + self.radius.setValue(0.1) + grid.addWidget(self.radius, 0, 1) # min no. locs - grid.addWidget(QtWidgets.QLabel("Min. no. locs:"), 2, 0) + grid.addWidget(QtWidgets.QLabel("Min. no. locs:"), 1, 0) self.min_locs = QtWidgets.QSpinBox() self.min_locs.setRange(1, int(1e6)) self.min_locs.setValue(10) - grid.addWidget(self.min_locs, 2, 1) + grid.addWidget(self.min_locs, 1, 1) # save cluster centers self.save_centers = QtWidgets.QCheckBox("Save cluster centers") self.save_centers.setChecked(False) - grid.addWidget(self.save_centers, 3, 0, 1, 2) + grid.addWidget(self.save_centers, 2, 0, 1, 2) # perform basic frame analysis self.frame_analysis = QtWidgets.QCheckBox( "Perform basic frame analysis" ) self.frame_analysis.setChecked(True) - grid.addWidget(self.frame_analysis, 4, 0, 1, 2) + grid.addWidget(self.frame_analysis, 3, 0, 1, 2) vbox.addLayout(grid) hbox = QtWidgets.QHBoxLayout() @@ -1978,31 +2050,33 @@ def __init__(self, window): def getParams(parent=None): """ Creates the dialog and returns the requested values for - SMLM clusterer (3D). + SMLM clusterer (2D). """ - dialog = SMLMDialog3D(parent) + dialog = SMLMDialog2D(parent) result = dialog.exec_() return ( - dialog.radius_xy.value(), - dialog.radius_z.value(), + dialog.radius.value(), dialog.min_locs.value(), dialog.save_centers.isChecked(), dialog.frame_analysis.isChecked(), result == QtWidgets.QDialog.Accepted, - ) + ) -class SMLMDialog2D(QtWidgets.QDialog): +class SMLMDialog3D(QtWidgets.QDialog): """ - A class to obtain inputs for SMLM clusterer (2D). + A class to obtain inputs for SMLM clusterer (3D). ... Attributes ---------- - radius : QDoubleSpinBox + radius_xy : QDoubleSpinBox contains clustering radius in x and y directions + radius_z : QDoubleSpinBox + contains clustering radius in z direction + (typically =2*radius_xy) min_locs : QSpinBox contains minimum number of locs in cluster @@ -2012,36 +2086,43 @@ class SMLMDialog2D(QtWidgets.QDialog): Creates the dialog and returns the requested values for clustering """ - + def __init__(self, window): super().__init__(window) self.window = window - self.setWindowTitle("Enter parameters (2D)") + self.setWindowTitle("Enter parameters (3D)") vbox = QtWidgets.QVBoxLayout(self) grid = QtWidgets.QGridLayout() - # clustering radius - grid.addWidget(QtWidgets.QLabel("Cluster radius (pixels):"), 0, 0) - self.radius = QtWidgets.QDoubleSpinBox() - self.radius.setRange(0.0001, 1e3) - self.radius.setDecimals(4) - self.radius.setValue(0.1) - grid.addWidget(self.radius, 0, 1) + # radius xy + grid.addWidget(QtWidgets.QLabel("Cluster radius xy (pixels):"), 0, 0) + self.radius_xy = QtWidgets.QDoubleSpinBox() + self.radius_xy.setRange(0.0001, 1e3) + self.radius_xy.setDecimals(4) + self.radius_xy.setValue(0.1) + grid.addWidget(self.radius_xy, 0, 1) + # radius z + grid.addWidget(QtWidgets.QLabel("Cluster radius z (pixels):"), 1, 0) + self.radius_z = QtWidgets.QDoubleSpinBox() + self.radius_z.setRange(0, 1e3) + self.radius_z.setDecimals(4) + self.radius_z.setValue(0.25) + grid.addWidget(self.radius_z, 1, 1) # min no. locs - grid.addWidget(QtWidgets.QLabel("Min. no. locs:"), 1, 0) + grid.addWidget(QtWidgets.QLabel("Min. no. locs:"), 2, 0) self.min_locs = QtWidgets.QSpinBox() self.min_locs.setRange(1, int(1e6)) self.min_locs.setValue(10) - grid.addWidget(self.min_locs, 1, 1) + grid.addWidget(self.min_locs, 2, 1) # save cluster centers self.save_centers = QtWidgets.QCheckBox("Save cluster centers") self.save_centers.setChecked(False) - grid.addWidget(self.save_centers, 2, 0, 1, 2) + grid.addWidget(self.save_centers, 3, 0, 1, 2) # perform basic frame analysis self.frame_analysis = QtWidgets.QCheckBox( "Perform basic frame analysis" ) self.frame_analysis.setChecked(True) - grid.addWidget(self.frame_analysis, 3, 0, 1, 2) + grid.addWidget(self.frame_analysis, 4, 0, 1, 2) vbox.addLayout(grid) hbox = QtWidgets.QHBoxLayout() @@ -2060,18 +2141,19 @@ def __init__(self, window): def getParams(parent=None): """ Creates the dialog and returns the requested values for - SMLM clusterer (2D). + SMLM clusterer (3D). """ - dialog = SMLMDialog2D(parent) + dialog = SMLMDialog3D(parent) result = dialog.exec_() return ( - dialog.radius.value(), + dialog.radius_xy.value(), + dialog.radius_z.value(), dialog.min_locs.value(), dialog.save_centers.isChecked(), dialog.frame_analysis.isChecked(), result == QtWidgets.QDialog.Accepted, - ) + ) class TestClustererDialog(QtWidgets.QDialog): @@ -5525,18 +5607,16 @@ class View(QtWidgets.QLabel): Called on pressing right arrow; moves FOV to_up() Called on pressing up arrow; moves FOV - undrift() - Undrifts with RCC - undrift_from_picked - Gets channel for undrifting from picked locs - _undrift_from_picked - Undrifts based on picked locs in a given channel - _undrift_from_picked_coordinate + undrift_aim() + Undrifts with AIM. + undrift_from_picked() + Undrifts from picked localizations. + _undrift_from_picked_coordinate() Calculates drift in a given coordinate - undrift_from_picked2d - Gets channel for undrifting from picked locs in 2D - _undrift_from_picked2d - Undrifts in x and y based on picked locs in a given channel + undrift_from_picked2d() + Undrifts x and y coordinates from picked localizations. + undrift_rcc() + Undrifts with RCC undo_drift Gets channel for undoing drift _undo_drift @@ -7063,7 +7143,7 @@ def get_channel(self, title="Choose a channel"): else: return None - def save_channel(self, title="Choose a channel"): + def save_channel(self, title="Choose a channel to save localizations"): """ Opens an input dialog to ask which channel to save. There is an option to save all channels separetely or merge @@ -7085,7 +7165,7 @@ def save_channel(self, title="Choose a channel"): pathlist.append("Combine all channels") index, ok = QtWidgets.QInputDialog.getItem( self, - "Save localizations", + title, "Channel:", pathlist, editable=False, @@ -7118,7 +7198,7 @@ def get_channel_all_seq(self, title="Choose a channel"): pathlist.append("Apply to all sequentially") index, ok = QtWidgets.QInputDialog.getItem( self, - "Save localizations", + title, "Channel:", pathlist, editable=False, @@ -9648,15 +9728,48 @@ def show_drift(self): self.plot_window.show() + def undrift_aim(self): + """Undrifts with Adaptive Intersection Maximization (AIM). + + See Ma H., et al. Science Advances. 2024.""" + + channel = self.get_channel("Undrift by AIM") + if channel is not None: + locs = self.all_locs[channel] + info = self.infos[channel] + pixelsize = self.window.display_settings_dlg.pixelsize.value() + + # get parameters for AIM + params, ok = AIMDialog.getParams() + params["intersect_d"] = params["intersect_d"] / pixelsize + params["roi_r"] = params["roi_r"] / pixelsize + if ok: + n_frames = info[0]["Frames"] + n_segments = int(np.floor(n_frames / params["segmentation"])) + progress = lib.ProgressDialog( + "Undrifting by AIM (1/2)", 0, n_segments, self.window + ) + locs, new_info, drift = aim.aim( + locs, info, **params, progress=progress + ) + # sanity check and assign attributes + locs = lib.ensure_sanity(locs, info) + self.all_locs[channel] = locs + self.locs[channel] = copy.copy(locs) + self.infos[channel] = new_info + self.index_blocks[channel] = None + self.add_drift(channel, drift) + self.update_scene() + self.show_drift() - def undrift(self): + def undrift_rcc(self): """ Undrifts with RCC. See Wang Y., et al. Optics Express. 2014 """ - channel = self.get_channel("Undrift") + channel = self.get_channel("Undrift by RCC") if channel is not None: info = self.infos[channel] n_frames = info[0]["Frames"] @@ -9721,23 +9834,78 @@ def undrift(self): @check_picks def undrift_from_picked(self): - """ Gets channel for undrifting from picked locs. """ + """Undrifts based on picked locs in a given channel.""" channel = self.get_channel("Undrift from picked") if channel is not None: - self._undrift_from_picked(channel) + picked_locs = self.picked_locs(channel) + status = lib.StatusDialog("Calculating drift...", self) + + drift_x = self._undrift_from_picked_coordinate( + channel, picked_locs, "x" + ) # find drift in x + drift_y = self._undrift_from_picked_coordinate( + channel, picked_locs, "y" + ) # find drift in y + + # Apply drift + self.all_locs[channel].x -= drift_x[self.all_locs[channel].frame] + self.all_locs[channel].y -= drift_y[self.all_locs[channel].frame] + self.locs[channel].x -= drift_x[self.locs[channel].frame] + self.locs[channel].y -= drift_y[self.locs[channel].frame] + + # A rec array to store the applied drift + drift = (drift_x, drift_y) + drift = np.rec.array(drift, dtype=[("x", "f"), ("y", "f")]) + + # If z coordinate exists, also apply drift there + if all([hasattr(_, "z") for _ in picked_locs]): + drift_z = self._undrift_from_picked_coordinate( + channel, picked_locs, "z" + ) + self.all_locs[channel].z -= drift_z[self.all_locs[channel].frame] + self.locs[channel].z -= drift_z[self.locs[channel].frame] + drift = lib.append_to_rec(drift, drift_z, "z") + + # Cleanup + self.index_blocks[channel] = None + self.add_drift(channel, drift) + status.close() + self.update_scene() @check_picks def undrift_from_picked2d(self): - """ - Gets channel for undrifting from picked locs in 2D. - Available when 3D data is loaded. - """ + """Undrifts in x and y based on picked locs in a given channel. + Available when 3D data is loaded.""" channel = self.get_channel("Undrift from picked") if channel is not None: - self._undrift_from_picked2d(channel) + picked_locs = self.picked_locs(channel) + status = lib.StatusDialog("Calculating drift...", self) + + drift_x = self._undrift_from_picked_coordinate( + channel, picked_locs, "x" + ) + drift_y = self._undrift_from_picked_coordinate( + channel, picked_locs, "y" + ) + # Apply drift + self.all_locs[channel].x -= drift_x[self.all_locs[channel].frame] + self.all_locs[channel].y -= drift_y[self.all_locs[channel].frame] + self.locs[channel].x -= drift_x[self.locs[channel].frame] + self.locs[channel].y -= drift_y[self.locs[channel].frame] + + # A rec array to store the applied drift + drift = (drift_x, drift_y) + drift = np.rec.array(drift, dtype=[("x", "f"), ("y", "f")]) + + # Cleanup + self.index_blocks[channel] = None + self.add_drift(channel, drift) + status.close() + self.update_scene() + def _undrift_from_picked_coordinate( self, channel, picked_locs, coordinate ): @@ -9795,87 +9963,6 @@ def nan_helper(y): return drift_mean - def _undrift_from_picked(self, channel): - """ - Undrifts based on picked locs in a given channel. - - Parameters - ---------- - channel : int - Channel to be undrifted - """ - - picked_locs = self.picked_locs(channel) - status = lib.StatusDialog("Calculating drift...", self) - - drift_x = self._undrift_from_picked_coordinate( - channel, picked_locs, "x" - ) # find drift in x - drift_y = self._undrift_from_picked_coordinate( - channel, picked_locs, "y" - ) # find drift in y - - # Apply drift - self.all_locs[channel].x -= drift_x[self.all_locs[channel].frame] - self.all_locs[channel].y -= drift_y[self.all_locs[channel].frame] - self.locs[channel].x -= drift_x[self.locs[channel].frame] - self.locs[channel].y -= drift_y[self.locs[channel].frame] - - # A rec array to store the applied drift - drift = (drift_x, drift_y) - drift = np.rec.array(drift, dtype=[("x", "f"), ("y", "f")]) - - # If z coordinate exists, also apply drift there - if all([hasattr(_, "z") for _ in picked_locs]): - drift_z = self._undrift_from_picked_coordinate( - channel, picked_locs, "z" - ) - self.all_locs[channel].z -= drift_z[self.all_locs[channel].frame] - self.locs[channel].z -= drift_z[self.locs[channel].frame] - drift = lib.append_to_rec(drift, drift_z, "z") - - # Cleanup - self.index_blocks[channel] = None - self.add_drift(channel, drift) - status.close() - self.update_scene() - - def _undrift_from_picked2d(self, channel): - """ - Undrifts in x and y based on picked locs in a given channel. - - Parameters - ---------- - channel : int - Channel to be undrifted - """ - - picked_locs = self.picked_locs(channel) - status = lib.StatusDialog("Calculating drift...", self) - - drift_x = self._undrift_from_picked_coordinate( - channel, picked_locs, "x" - ) - drift_y = self._undrift_from_picked_coordinate( - channel, picked_locs, "y" - ) - - # Apply drift - self.all_locs[channel].x -= drift_x[self.all_locs[channel].frame] - self.all_locs[channel].y -= drift_y[self.all_locs[channel].frame] - self.locs[channel].x -= drift_x[self.locs[channel].frame] - self.locs[channel].y -= drift_y[self.locs[channel].frame] - - # A rec array to store the applied drift - drift = (drift_x, drift_y) - drift = np.rec.array(drift, dtype=[("x", "f"), ("y", "f")]) - - # Cleanup - self.index_blocks[channel] = None - self.add_drift(channel, drift) - status.close() - self.update_scene() - def undo_drift(self): """ Gets channel for undoing drift. """ @@ -10900,9 +10987,9 @@ def initUI(self, plugins_loaded): # menu bar - Postprocess postprocess_menu = self.menu_bar.addMenu("Postprocess") - undrift_action = postprocess_menu.addAction("Undrift by RCC") - undrift_action.setShortcut("Ctrl+U") - undrift_action.triggered.connect(self.view.undrift) + undrift_aim_action = postprocess_menu.addAction("Undrift by AIM") + undrift_aim_action.setShortcut("Ctrl+U") + undrift_aim_action.triggered.connect(self.view.undrift_aim) undrift_from_picked_action = postprocess_menu.addAction( "Undrift from picked" ) @@ -10916,6 +11003,8 @@ def initUI(self, plugins_loaded): undrift_from_picked2d_action.triggered.connect( self.view.undrift_from_picked2d ) + undrift_action = postprocess_menu.addAction("Undrift by RCC") + undrift_action.triggered.connect(self.view.undrift_rcc) drift_action = postprocess_menu.addAction("Undo drift") drift_action.triggered.connect(self.view.undo_drift) drift_action = postprocess_menu.addAction("Show drift") diff --git a/picasso/imageprocess.py b/picasso/imageprocess.py index 892488b5..6b63aefd 100644 --- a/picasso/imageprocess.py +++ b/picasso/imageprocess.py @@ -1,5 +1,5 @@ """ - picasso/imageprocess + picasso.imageprocess ~~~~~~~~~~~~~~~~~~~~ Image processing functions diff --git a/picasso/lib.py b/picasso/lib.py index bae337bc..2e0e6be1 100644 --- a/picasso/lib.py +++ b/picasso/lib.py @@ -57,6 +57,13 @@ def set_value(self, value): def closeEvent(self, event): _dialogs.remove(self) + def zero_progress(self, description=None): + """Sets progress dialog to zero and changes title if given.""" + + if description: + self.setLabelText(description) + self.set_value(0) + class StatusDialog(QtWidgets.QDialog): """StatusDialog displays the description string in a dialog.""" diff --git a/readme.rst b/readme.rst index 1374679a..b17a676e 100644 --- a/readme.rst +++ b/readme.rst @@ -24,13 +24,9 @@ A collection of tools for painting super-resolution images. The Picasso software A comprehensive documentation can be found here: `Read the Docs `__. -Picasso 0.6.1 +Picasso 0.7.0 ------------- -In previous version, the rotation window (3D Render) showed an incorrect length of the scalebar. This has been fixed. - -Picasso 0.6.0 -------------- -RESI dialog added to Picasso Render, allowing for substantial boost in spatial resolution (*Reinhardt, et al., Nature, 2023.* DOI: 10.1038/s41586-023-05925-9). +Adaptive Intersection Maximization (AIM, `doi: 10.1126/sciadv.adm7765 `_) implemented in Picasso. The implementation may not be as fast as in the original version due to Python constraints. Previous versions ----------------- From cd4b8e7ac536b6396c11a0f101a30da4ab80d513 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 5 Jul 2024 14:49:56 +0200 Subject: [PATCH 04/16] CLEAN --- picasso/aim.py | 2 -- picasso/gui/render.py | 6 ------ 2 files changed, 8 deletions(-) diff --git a/picasso/aim.py b/picasso/aim.py index 650969ef..dbbe2074 100644 --- a/picasso/aim.py +++ b/picasso/aim.py @@ -19,8 +19,6 @@ _InterpolatedUnivariateSpline from tqdm import tqdm as _tqdm -#TODO: add progress dialog, also tqdm - def intersect1d(a, b): """Slightly faster implementation of _np.intersect1d without diff --git a/picasso/gui/render.py b/picasso/gui/render.py index 4600c954..950a3a02 100644 --- a/picasso/gui/render.py +++ b/picasso/gui/render.py @@ -9795,7 +9795,6 @@ def undrift_rcc(self): "Correlating image pairs", 0, n_pairs, self ) try: - start_time = time.time() # find drift and apply it to locs drift, _ = postprocess.undrift( locs, @@ -9805,11 +9804,6 @@ def undrift_rcc(self): seg_progress.set_value, rcc_progress.set_value, ) - finish_time = time.time() - print( - "RCC drift estimate running time [seconds]: ", - np.round(finish_time-start_time, 1) - ) # sanity check and assign attributes locs = lib.ensure_sanity(locs, info) self.all_locs[channel] = locs From 81242938ed7505c988a643d3a55db972eb2b704a Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 5 Jul 2024 14:54:30 +0200 Subject: [PATCH 05/16] rotation window from rectangular and polygonal pick FIX --- changelog.rst | 1 + picasso/gui/rotation.py | 9 ++++----- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/changelog.rst b/changelog.rst index 7a2a3d7d..f8f6824f 100644 --- a/changelog.rst +++ b/changelog.rst @@ -7,6 +7,7 @@ Last change: 28-JUN-2024 MTS ----- - Adaptive Intersection Maximization (AIM, doi: 10.1038/s41592-022-01307-0) implemented - CMD ``clusterfile`` fixed +- Picasso: Render 3D, rectangular and polygonal pick fixed 0.6.9 - 0.6.11 -------------- diff --git a/picasso/gui/rotation.py b/picasso/gui/rotation.py index b4a8a96a..4f61c818 100644 --- a/picasso/gui/rotation.py +++ b/picasso/gui/rotation.py @@ -22,8 +22,7 @@ from numpy.lib.recfunctions import stack_arrays -from .. import io, render -from ..lib import ProgressDialog +from .. import io, render, lib DEFAULT_OVERSAMPLING = 1.0 @@ -614,7 +613,7 @@ def build_animation(self): # render all frames and save in RAM video_writer = imageio.get_writer(name, fps=self.fps.value()) - progress = ProgressDialog( + progress = lib.ProgressDialog( "Rendering frames", 0, len(angx), self.window ) progress.set_value(0) @@ -1528,7 +1527,7 @@ def fit_in_view_rotated(self, get_viewport=False): elif self.pick_shape == "Rectangle": w = self.pick_size (xs, ys), (xe, ye) = self.pick - X, Y = self.window.window.view.get_pick_rectangle_corners( + X, Y = lib.get_pick_rectangle_corners( xs, ys, xe, ye, w ) x_min = min(X) @@ -1536,7 +1535,7 @@ def fit_in_view_rotated(self, get_viewport=False): y_min = min(Y) y_max = max(Y) elif self.pick_shape == "Polygon": - X, Y = self.window.window.view.get_pick_polygon_corners(self.pick) + X, Y = lib.get_pick_polygon_corners(self.pick) x_min = min(X) x_max = max(X) y_min = min(Y) From 2d810ccb0a80bdb02864945ff1f173b8ffb0820e Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 5 Jul 2024 15:22:54 +0200 Subject: [PATCH 06/16] Fit Z, fix bug with overshooting minimize_scalar to avoid gaps in calibration --- picasso/zfit.py | 6 +++++- readme.rst | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/picasso/zfit.py b/picasso/zfit.py index 082a4d2b..9bc39da5 100644 --- a/picasso/zfit.py +++ b/picasso/zfit.py @@ -220,7 +220,11 @@ def fit_z( sx = locs.sx sy = locs.sy for i in range(len(z)): - result = _minimize_scalar(_fit_z_target, args=(sx[i], sy[i], cx, cy)) + result = _minimize_scalar( + _fit_z_target, + bounds=[-1000, 1000], # to avoid potential gaps in the calibration curve, credits to Loek Andriessen + args=(sx[i], sy[i], cx, cy) + ) z[i] = result.x square_d_zcalib[i] = result.fun z *= magnification_factor diff --git a/readme.rst b/readme.rst index b17a676e..a4eec3b3 100644 --- a/readme.rst +++ b/readme.rst @@ -26,7 +26,7 @@ A comprehensive documentation can be found here: `Read the Docs `_) implemented in Picasso. The implementation may not be as fast as in the original version due to Python constraints. +Adaptive Intersection Maximization (AIM, `doi: 10.1126/sciadv.adm7765 `_) implemented in Picasso. Previous versions ----------------- From 58aa7626056d982430920b3970d61eba02a80e5e Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Fri, 5 Jul 2024 16:11:02 +0200 Subject: [PATCH 07/16] default MLE fitting should uses sigmaxy method (not sigma) --- picasso/localize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/picasso/localize.py b/picasso/localize.py index 89de078c..a5bffe8d 100755 --- a/picasso/localize.py +++ b/picasso/localize.py @@ -391,7 +391,7 @@ def fit_async( box, eps=0.001, max_it=100, - method="sigma", + method="sigmaxy", ): spots = get_spots(movie, identifications, box, camera_info) return _gaussmle.gaussmle_async(spots, eps, max_it, method=method) From 7a7cbefdba69f39d91ec80eb81064280377eac9f Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 17:53:48 +0200 Subject: [PATCH 08/16] picasso.localize.localize FIX --- changelog.rst | 5 ++++- picasso/localize.py | 10 +++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/changelog.rst b/changelog.rst index f8f6824f..9eb18fe7 100644 --- a/changelog.rst +++ b/changelog.rst @@ -1,13 +1,16 @@ Changelog ========= -Last change: 28-JUN-2024 MTS +Last change: 06-JUL-2024 MTS 0.7.0 ----- - Adaptive Intersection Maximization (AIM, doi: 10.1038/s41592-022-01307-0) implemented +- Z fitting improved by setting bounds on fitted z values to avoid NaNs - CMD ``clusterfile`` fixed - Picasso: Render 3D, rectangular and polygonal pick fixed +- picasso.localize.localize fixed +- default MLE fitting uses different sx and sy (CMD only) 0.6.9 - 0.6.11 -------------- diff --git a/picasso/localize.py b/picasso/localize.py index a5bffe8d..bc3c0b23 100755 --- a/picasso/localize.py +++ b/picasso/localize.py @@ -424,9 +424,13 @@ def locs_from_fits(identifications, theta, CRLBs, likelihoods, iterations, box): return locs -def localize(movie, info, parameters): - identifications = identify(movie, parameters) - return fit(movie, info, identifications, parameters["Box Size"]) +def localize(movie, camera_info, parameters): + identifications = identify( + movie, + parameters["Min. Net Gradient"], + parameters["Box Size"], + ) + return fit(movie, camera_info, identifications, parameters["Box Size"]) def check_nena(locs, info, callback=None): From 99c8b6e97e2fb959f27d80fb6d16f5718b158134 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 18:08:37 +0200 Subject: [PATCH 09/16] AIM documentation --- docs/render.rst | 29 ++++++++++++++++++++++------- picasso/gui/render.py | 2 +- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/docs/render.rst b/docs/render.rst index b5e182e3..7386c3eb 100644 --- a/docs/render.rst +++ b/docs/render.rst @@ -14,13 +14,16 @@ Opening Files Drift Correction ---------------- -Picasso offers two procedures to correct for drift: an RCC algorithm (option A), and use of specific structures in the image as drift markers (option B). Although option A does not require any additional sample preparation, option B depends on the presence of either fiducial markers or inherently clustered structures in the image. On the other hand, option B often supports more precise drift estimation and thus allows for higher image resolution. To achieve the highest possible resolution (ultra-resolution), we recommend consecutive applications of option A and multiple rounds of option B. The drift markers for option B can be features of the image itself (e.g., protein complexes or DNA origami) or intentionally included markers (e.g., DNA origami or gold nanoparticles). When using DNA origami as drift markers, the correction is typically applied in two rounds: first, with whole DNA origami structures as markers, and, second, using single DNA-PAINT binding sites as markers. In both cases, the precision of drift correction strongly depends on the number of selected drift markers. +Picasso offers three procedures to correct for drift: AIM (Ma, H., et al. Science Advances. 2024., option A), use of specific structures in the image as drift markers (option B) and an RCC algorithm (option C). AIM is precise, robust, quick, requires no user interaction or fiducial markers (although adding them will may improve performance). Although RCC does not require any additional sample preparation, option B depends on the presence of either fiducial markers or inherently clustered structures in the image. On the other hand, option B often supports more precise drift estimation and thus allows for higher image resolution. To achieve the highest possible resolution (ultra-resolution), we recommend AIM or consecutive applications of option C and multiple rounds of option B. The drift markers for option B can be features of the image itself (e.g., protein complexes or DNA origami) or intentionally included markers (e.g., DNA origami or gold nanoparticles). When using DNA origami as drift markers, the correction is typically applied in two rounds: first, with whole DNA origami structures as markers, and, second, using single DNA-PAINT binding sites as markers. In both cases, the precision of drift correction strongly depends on the number of selected drift markers. -Redundant cross-correlation drift correction -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Adaptive Intersection Maximization (AIM) drift correction +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -1. In ``Picasso: Render``, select ``Postprocess > Undrift by RCC``. -2. A dialog will appear asking for the segmentation parameter. Although the default value, 1,000 frames, is a sensible choice for most movies, it might be necessary to adjust the segmentation parameter of the algorithm, depending on the total number of frames in the movie and the number of localizations per frame. A smaller segment size results in better temporal drift resolution but requires a movie with more localizations per frame. +1. In ``Picasso: Render``, select ``Postprocess > Undrift by AIM``. +2. The dialog asks the user to select: + a. ``Segmentation`` - the number of frames per interval to calculate the drift. The lower the value, the better the temporal resolution of the drift correction, but the higher the computational cost. + b. ``Intersection distance (nm)`` - the maximum distance between two localizations in two consecutive temporal segments to be considered the same molecule. This parameter is robust, 3*NeNA for optimal result is recommended. + c. ``Max. drift in segment (nm)`` - the maximum expected drift between two consecutive temporal segments. If the drift is larger, the algorithm will likely diverge. Setting the parameter up to ``3 * intersection_distance`` will result in fast computation. 3. After the algorithm finishes, the estimated drift will be displayed in a pop-up window, and the display will show the drift-corrected image. Marker-based drift correction @@ -31,6 +34,14 @@ Marker-based drift correction 3. Select ``Postprocess > Undrift from picked`` to compute and apply the drift correction. 4. (Optional) Save the drift-corrected localizations by selecting ``File > Save localizations``. +Redundant cross-correlation drift correction +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +1. In ``Picasso: Render``, select ``Postprocess > Undrift by RCC``. +2. A dialog will appear asking for the segmentation parameter. Although the default value, 1,000 frames, is a sensible choice for most movies, it might be necessary to adjust the segmentation parameter of the algorithm, depending on the total number of frames in the movie and the number of localizations per frame. A smaller segment size results in better temporal drift resolution but requires a movie with more localizations per frame. +3. After the algorithm finishes, the estimated drift will be displayed in a pop-up window, and the display will show the drift-corrected image. + + Picking of regions of interest ------------------------------ @@ -340,9 +351,9 @@ Allows the user to display only a fraction of localizations to speed up renderin Postprocess ~~~~~~~~~~~ -Undrift by RCC +Undrift by AIM ^^^^^^^^^^^^^^ -Performs drift correction by redundant cross-correlation. +Performs drift correction using the AIM algorithm (Ma, H., et al. Science Advances. 2024). Undrift from picked (3D) ^^^^^^^^^^^^^^^^^^^^^^^^ @@ -352,6 +363,10 @@ Undrift from picked (2D) ^^^^^^^^^^^^^^^^^^^^^^^^ Performs drift correction using the picked localizations as fiducials. Does not perform drift correction in z even if dataset has 3D information. +Undrift by RCC +^^^^^^^^^^^^^^ +Performs drift correction by redundant cross-correlation. + Undo drift (2D) ^^^^^^^^^^^^^^^ Undo previous drift correction (only 2D part). Can be pressed again to redo. diff --git a/picasso/gui/render.py b/picasso/gui/render.py index 950a3a02..db8af49f 100644 --- a/picasso/gui/render.py +++ b/picasso/gui/render.py @@ -1739,7 +1739,7 @@ def __init__(self, window): self.intersect_d.setSingleStep(1) grid.addWidget(self.intersect_d, 1, 1) grid.addWidget( - QtWidgets.QLabel("Max. drift in segmentation (nm):"), 2, 0 + QtWidgets.QLabel("Max. drift in segment (nm):"), 2, 0 ) self.max_drift = QtWidgets.QDoubleSpinBox() self.max_drift.setRange(0.1, 1e6) From 6dfd1b284c84d1e6f398db5c6bf4b89c2d1820bc Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 18:23:03 +0200 Subject: [PATCH 10/16] citation advice for external software (NeNA, AIM, etc) --- readme.rst | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/readme.rst b/readme.rst index a4eec3b3..29c4e863 100644 --- a/readme.rst +++ b/readme.rst @@ -122,7 +122,7 @@ Contributions & Copyright | Contributors: Joerg Schnitzbauer, Maximilian Strauss, Rafal Kowalewski, Adrian Przybylski, Andrey Aristov, Hiroshi Sasaki, Alexander Auer, Johanna Rahm | Copyright (c) 2015-2019 Jungmann Lab, Max Planck Institute of Biochemistry | Copyright (c) 2020-2021 Maximilian Strauss -| Copyright (c) 2022-2023 Rafal Kowalewski +| Copyright (c) 2022-2024 Rafal Kowalewski Citing Picasso -------------- @@ -133,6 +133,13 @@ If you use picasso in your research, please cite our Nature Protocols publicatio | Super-Resolution Microscopy with DNA-PAINT | Nature Protocols (2017). 12: 1198-1228 DOI: `https://doi.org/10.1038/nprot.2017.024 `__ +If you use some of the functionalities provided by Picasso, please also cite the respective publications: + +- Nearest Neighbor based Analysis (NeNA) for experimental localization precision. DOI: `https://doi.org/10.1007/s00418-014-1192-3 `__ +- Theoretical localization precision (Gauss LQ and MLE). DOI: `https://doi.org/10.1038/nmeth.1447 `__ +- MLE fitting. DOI: `https://doi.org/10.1038/nmeth.1449 `__ +- AIM undrifting. DOI: `10.1126/sciadv.adm776 `__ + Credits ------- From b9b8d2e9876f6b4e8face4fde85494b9cf1ce4e5 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 18:29:02 +0200 Subject: [PATCH 11/16] CLEAN --- readme.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/readme.rst b/readme.rst index 29c4e863..79a45977 100644 --- a/readme.rst +++ b/readme.rst @@ -133,6 +133,7 @@ If you use picasso in your research, please cite our Nature Protocols publicatio | Super-Resolution Microscopy with DNA-PAINT | Nature Protocols (2017). 12: 1198-1228 DOI: `https://doi.org/10.1038/nprot.2017.024 `__ + If you use some of the functionalities provided by Picasso, please also cite the respective publications: - Nearest Neighbor based Analysis (NeNA) for experimental localization precision. DOI: `https://doi.org/10.1007/s00418-014-1192-3 `__ From 6b2828b4ec5d934bf2cdef892af2dfe648eb114b Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 18:30:27 +0200 Subject: [PATCH 12/16] CLEAN --- readme.rst | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/readme.rst b/readme.rst index 79a45977..eee8a860 100644 --- a/readme.rst +++ b/readme.rst @@ -133,8 +133,7 @@ If you use picasso in your research, please cite our Nature Protocols publicatio | Super-Resolution Microscopy with DNA-PAINT | Nature Protocols (2017). 12: 1198-1228 DOI: `https://doi.org/10.1038/nprot.2017.024 `__ - -If you use some of the functionalities provided by Picasso, please also cite the respective publications: +| If you use some of the functionalities provided by Picasso, please also cite the respective publications: - Nearest Neighbor based Analysis (NeNA) for experimental localization precision. DOI: `https://doi.org/10.1007/s00418-014-1192-3 `__ - Theoretical localization precision (Gauss LQ and MLE). DOI: `https://doi.org/10.1038/nmeth.1447 `__ From 1dc33f5445963d830e41e782f1fbe07fb5294bbd Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Sat, 6 Jul 2024 18:32:08 +0200 Subject: [PATCH 13/16] CLEAN --- readme.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.rst b/readme.rst index eee8a860..ea96bea2 100644 --- a/readme.rst +++ b/readme.rst @@ -132,7 +132,7 @@ If you use picasso in your research, please cite our Nature Protocols publicatio | J. Schnitzbauer*, M.T. Strauss*, T. Schlichthaerle, F. Schueder, R. Jungmann | Super-Resolution Microscopy with DNA-PAINT | Nature Protocols (2017). 12: 1198-1228 DOI: `https://doi.org/10.1038/nprot.2017.024 `__ - +| | If you use some of the functionalities provided by Picasso, please also cite the respective publications: - Nearest Neighbor based Analysis (NeNA) for experimental localization precision. DOI: `https://doi.org/10.1007/s00418-014-1192-3 `__ From 82ff81cd94920ceae7e19d6546ea72caa6f7366d Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Mon, 8 Jul 2024 16:06:01 +0200 Subject: [PATCH 14/16] fix the segmentation boundaries in AIM; 2nd round undrifts the first interval too --- picasso/aim.py | 96 +++++++++++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 41 deletions(-) diff --git a/picasso/aim.py b/picasso/aim.py index dbbe2074..d7ccf61f 100644 --- a/picasso/aim.py +++ b/picasso/aim.py @@ -353,7 +353,8 @@ def get_fft_peak_z(roi_cc, roi_size): def intersection_max( x, y, ref_x, ref_y, - frame, segmentation, intersect_d, roi_r, width, progress=None, + frame, seg_bounds, intersect_d, roi_r, width, + aim_round=1, progress=None, ): """Maximize intersection (undrift) for 2D localizations. @@ -369,8 +370,9 @@ def intersection_max( y-coordinates of the reference localizations. frame : _np.array Frame indices of localizations. - segmentation : int - Time interval for drift tracking, unit: frames. + seg_bounds : np.array + Frame indices of the segmentation bounds. Defines temporal + intervals used to estimate drift. intersect_d : float Intersect distance in camera pixels. roi_r : float @@ -378,6 +380,11 @@ def intersection_max( higher than the maximum expected drift within one segment. width : int Width of the camera image in camera pixels. + aim_round : {1, 2} + Round of AIM algorithm. The first round uses the first interval + as reference, the second round uses the entire dataset as + reference. The impact is that in the second round, the first + interval is also undrifted. progress : picasso.lib.ProgressDialog (default=None) Progress dialog. If None, progress is displayed with tqdm. @@ -393,8 +400,10 @@ def intersection_max( Drift in y-direction. """ + assert aim_round in [1, 2], "aim_round must be 1 or 2." + # number of segments - n_segments = int(_np.ceil(frame.max() / segmentation)) + n_segments = len(seg_bounds) - 1 rel_drift_x = 0 # adaptive drift (updated at each interval) rel_drift_y = 0 @@ -411,7 +420,7 @@ def intersection_max( for i, shift_x in enumerate(steps): for j, shift_y in enumerate(steps): shifts_xy[i, j] = shift_x + shift_y * width_units - shifts_xy = shifts_xy.reshape(box**2) + shifts_xy = shifts_xy.reshape(box ** 2) # convert reference to a 1D array in units of intersect_d and find # unique values and counts @@ -421,16 +430,21 @@ def intersection_max( l0_coords, l0_counts = _np.unique(l0, return_counts=True) # initialize progress such that if GUI is used, tqdm is omitted + start_idx = 1 if aim_round == 1 else 0 if progress is not None: - iterator = range(1, n_segments) + iterator = range(start_idx, n_segments) else: - iterator = _tqdm(n_segments, desc="Undrifting z", unit="segment") + iterator = _tqdm( + range(start_idx, n_segments), + desc=f"Undrifting ({aim_round}/2)", + unit="segment", + ) # run across each segment for s in iterator: # get the target localizations within the current segment - min_frame_idx = frame > s * segmentation - max_frame_idx = frame <= (s + 1) * segmentation + min_frame_idx = frame > seg_bounds[s] + max_frame_idx = frame <= seg_bounds[s+1] x1 = x[min_frame_idx & max_frame_idx] y1 = y[min_frame_idx & max_frame_idx] @@ -462,12 +476,10 @@ def intersection_max( iterator.update(s - iterator.n) # interpolate the drifts (cubic spline) for all frames - n_frames = n_segments * segmentation - drift_track_points = _np.linspace(0, n_frames, n_segments+1) - t = (drift_track_points[1:] + drift_track_points[:-1]) / 2 + t = (seg_bounds[1:] + seg_bounds[:-1]) / 2 drift_x_pol = _InterpolatedUnivariateSpline(t, drift_x, k=3) drift_y_pol = _InterpolatedUnivariateSpline(t, drift_y, k=3) - t_inter = _np.arange(n_frames) + 1 + t_inter = _np.arange(seg_bounds[-1]) + 1 drift_x = drift_x_pol(t_inter) drift_y = drift_y_pol(t_inter) @@ -480,8 +492,8 @@ def intersection_max( def intersection_max_z( x, y, z, ref_x, ref_y, ref_z, - frame, segmentation, intersect_d, roi_r, width, height, pixelsize, - progress=None, + frame, seg_bounds, intersect_d, roi_r, width, height, pixelsize, + aim_round=1, progress=None, ): """Maximize intersection (undrift) for 3D localizations. Assumes that x and y coordinates were already undrifted. x and y are in @@ -494,7 +506,7 @@ def intersection_max_z( ref_z = ref_z.copy() / pixelsize #TODO: remember to convert back to nm, also for drift_z # number of segments - n_segments = int(_np.ceil(frame.max() / segmentation)) + n_segments = len(seg_bounds) - 1 rel_drift_z = 0 # adaptive drift (updated at each interval) # drift in z @@ -521,16 +533,21 @@ def intersection_max_z( l0_coords, l0_counts = _np.unique(l0, return_counts=True) # initialize progress such that if GUI is used, tqdm is omitted + start_idx = 1 if aim_round == 1 else 0 if progress is not None: - iterator = range(1, n_segments) + iterator = range(start_idx, n_segments) else: - iterator = _tqdm(n_segments, desc="Undrifting z", unit="segment") + iterator = _tqdm( + range(start_idx, n_segments), + desc=f"Undrifting z ({aim_round}/2)", + unit="segment", + ) # run across each segment for s in iterator: # get the target localizations within the current segment - min_frame_idx = frame > s * segmentation - max_frame_idx = frame <= (s + 1) * segmentation + min_frame_idx = frame > seg_bounds[s] + max_frame_idx = frame <= seg_bounds[s+1] x1 = x[min_frame_idx & max_frame_idx] y1 = y[min_frame_idx & max_frame_idx] z1 = z[min_frame_idx & max_frame_idx] @@ -561,11 +578,9 @@ def intersection_max_z( # interpolate the drifts (cubic spline) for all frames - n_frames = n_segments * segmentation - drift_track_points = _np.linspace(0, n_frames, n_segments+1) - t = (drift_track_points[1:] + drift_track_points[:-1]) / 2 + t = (seg_bounds[1:] + seg_bounds[:-1]) / 2 drift_z_pol = _InterpolatedUnivariateSpline(t, drift_z, k=3) - t_inter = _np.arange(n_frames) + 1 + t_inter = _np.arange(seg_bounds[-1]) + 1 drift_z = drift_z_pol(t_inter) # undrift the localizations @@ -616,15 +631,14 @@ def aim( width = info[0]["Width"] height = info[0]["Height"] pixelsize = info[1]["Pixelsize"] + n_frames = info[0]["Frames"] # frames should start at 1 frame = locs["frame"] + 1 - n_frames = _np.max(frame) - - # group frames into track intervals and correct for rounding errors - n_segments = _np.floor(n_frames / segmentation).astype(int) - n_frames = n_segments * segmentation - frame[frame > n_frames] = n_frames # cap frame number + # find the segmentation bounds (temporal intervals) + seg_bounds = _np.concatenate(( + _np.arange(0, n_frames, segmentation), [n_frames] + )) # get the reference localizations (first interval) ref_x = locs["x"][frame <= segmentation] @@ -634,25 +648,25 @@ def aim( # the first run is with the first interval as reference x_pdc, y_pdc, drift_x1, drift_y1 = intersection_max( locs.x, locs.y, ref_x, ref_y, - frame, segmentation, intersect_d, roi_r, width, - progress=progress, + frame, seg_bounds, intersect_d, roi_r, width, + aim_round=1, progress=progress, ) # the second run is with the entire dataset as reference if progress is not None: progress.zero_progress(description="Undrifting by AIM (2/2)") x_pdc, y_pdc, drift_x2, drift_y2 = intersection_max( x_pdc, y_pdc, x_pdc, y_pdc, - frame, segmentation, intersect_d, roi_r, width, - progress=progress, + frame, seg_bounds, intersect_d, roi_r, width, + aim_round=2, progress=progress, ) # add the drifts together from the two rounds drift_x = drift_x1 + drift_x2 drift_y = drift_y1 + drift_y2 - # shift the drifts by the mean value (like in Picasso) - drift_x -= _np.mean(drift_x) - drift_y -= _np.mean(drift_y) + # # shift the drifts by the mean value (like in Picasso) + # drift_x -= _np.mean(drift_x) + # drift_y -= _np.mean(drift_y) # combine to Picasso format drift = _np.rec.array((drift_x, drift_y), dtype=[("x", "f"), ("y", "f")]) @@ -666,15 +680,15 @@ def aim( ref_z = locs.z[frame <= segmentation] z_pdc, drift_z1 = intersection_max_z( x_pdc, y_pdc, locs.z, ref_x, ref_y, ref_z, - frame, segmentation, intersect_d, roi_r, width, height, pixelsize, - progress=progress, + frame, seg_bounds, intersect_d, roi_r, width, height, pixelsize, + aim_round=1, progress=progress, ) if progress is not None: progress.zero_progress(description="Undrifting z (2/2)") z_pdc, drift_z2 = intersection_max_z( x_pdc, y_pdc, z_pdc, x_pdc, y_pdc, z_pdc, - frame, segmentation, intersect_d, roi_r, width, height, pixelsize, - progress=progress, + frame, seg_bounds, intersect_d, roi_r, width, height, pixelsize, + aim_round=2, progress=progress, ) drift_z = drift_z1 + drift_z2 drift_z -= _np.mean(drift_z) From e5afd2b8cf1a13fec66b2fc56ab8f2c3a87156e0 Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Mon, 8 Jul 2024 18:08:23 +0200 Subject: [PATCH 15/16] FIX --- picasso/aim.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/picasso/aim.py b/picasso/aim.py index d7ccf61f..5939e591 100644 --- a/picasso/aim.py +++ b/picasso/aim.py @@ -448,6 +448,12 @@ def intersection_max( x1 = x[min_frame_idx & max_frame_idx] y1 = y[min_frame_idx & max_frame_idx] + # skip if no reference localizations + if len(x1) == 0: + drift_x[s] = drift_x[s-1] + drift_y[s] = drift_y[s-1] + continue + # undrifting from the previous round x1 += rel_drift_x y1 += rel_drift_y @@ -552,6 +558,11 @@ def intersection_max_z( y1 = y[min_frame_idx & max_frame_idx] z1 = z[min_frame_idx & max_frame_idx] + # skip if no reference localizations + if len(x1) == 0: + drift_z[s] = drift_z[s-1] + continue + # undrifting from the previous round z1 += rel_drift_z @@ -664,9 +675,9 @@ def aim( drift_x = drift_x1 + drift_x2 drift_y = drift_y1 + drift_y2 - # # shift the drifts by the mean value (like in Picasso) - # drift_x -= _np.mean(drift_x) - # drift_y -= _np.mean(drift_y) + # # shift the drifts by the mean value + drift_x -= _np.mean(drift_x) + drift_y -= _np.mean(drift_y) # combine to Picasso format drift = _np.rec.array((drift_x, drift_y), dtype=[("x", "f"), ("y", "f")]) From 600a3d7e89b57554555af0d02e0c6708bb5d5f2f Mon Sep 17 00:00:00 2001 From: rafalkowalewski1 Date: Mon, 8 Jul 2024 19:01:00 +0200 Subject: [PATCH 16/16] =?UTF-8?q?Bump=20version:=200.6.11=20=E2=86=92=200.?= =?UTF-8?q?7.0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- distribution/picasso.iss | 4 ++-- docs/conf.py | 2 +- picasso/__init__.py | 2 +- picasso/__version__.py | 2 +- release/one_click_windows_gui/create_installer_windows.bat | 2 +- release/one_click_windows_gui/picasso_innoinstaller.iss | 4 ++-- setup.py | 2 +- 8 files changed, 10 insertions(+), 10 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index d17b37ca..8dba192c 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.6.11 +current_version = 0.7.0 commit = True tag = False parse = (?P\d+)\.(?P\d+)\.(?P\d+)(\-(?P[a-z]+)(?P\d+))? diff --git a/distribution/picasso.iss b/distribution/picasso.iss index ea0bb944..48643090 100644 --- a/distribution/picasso.iss +++ b/distribution/picasso.iss @@ -2,10 +2,10 @@ AppName=Picasso AppPublisher=Jungmann Lab, Max Planck Institute of Biochemistry -AppVersion=0.6.11 +AppVersion=0.7.0 DefaultDirName={commonpf}\Picasso DefaultGroupName=Picasso -OutputBaseFilename="Picasso-Windows-64bit-0.6.11" +OutputBaseFilename="Picasso-Windows-64bit-0.7.0" ArchitecturesAllowed=x64 ArchitecturesInstallIn64BitMode=x64 diff --git a/docs/conf.py b/docs/conf.py index 85037a37..ca5a3bb8 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -26,7 +26,7 @@ # The short X.Y version version = "" # The full version, including alpha/beta/rc tags -release = "0.6.11" +release = "0.7.0" # -- General configuration --------------------------------------------------- diff --git a/picasso/__init__.py b/picasso/__init__.py index ea75779a..b6273df9 100644 --- a/picasso/__init__.py +++ b/picasso/__init__.py @@ -8,7 +8,7 @@ import os.path as _ospath import yaml as _yaml -__version__ = "0.6.11" +__version__ = "0.7.0" _this_file = _ospath.abspath(__file__) _this_dir = _ospath.dirname(_this_file) diff --git a/picasso/__version__.py b/picasso/__version__.py index b9326a7c..c03d0cb3 100644 --- a/picasso/__version__.py +++ b/picasso/__version__.py @@ -1 +1 @@ -VERSION_NO = "0.6.11" +VERSION_NO = "0.7.0" diff --git a/release/one_click_windows_gui/create_installer_windows.bat b/release/one_click_windows_gui/create_installer_windows.bat index 9e33d299..08d0a8b2 100644 --- a/release/one_click_windows_gui/create_installer_windows.bat +++ b/release/one_click_windows_gui/create_installer_windows.bat @@ -11,7 +11,7 @@ call conda activate picasso_installer call python setup.py sdist bdist_wheel call cd release/one_click_windows_gui -call pip install "../../dist/picassosr-0.6.11-py3-none-any.whl" +call pip install "../../dist/picassosr-0.7.0-py3-none-any.whl" call pip install pyinstaller==5.7 call pyinstaller ../pyinstaller/picasso.spec -y --clean diff --git a/release/one_click_windows_gui/picasso_innoinstaller.iss b/release/one_click_windows_gui/picasso_innoinstaller.iss index e2b753f4..1165f1a1 100644 --- a/release/one_click_windows_gui/picasso_innoinstaller.iss +++ b/release/one_click_windows_gui/picasso_innoinstaller.iss @@ -1,10 +1,10 @@ [Setup] AppName=Picasso AppPublisher=Jungmann Lab, Max Planck Institute of Biochemistry -AppVersion=0.6.11 +AppVersion=0.7.0 DefaultDirName={commonpf}\Picasso DefaultGroupName=Picasso -OutputBaseFilename="Picasso-Windows-64bit-0.6.11" +OutputBaseFilename="Picasso-Windows-64bit-0.7.0" ArchitecturesAllowed=x64 ArchitecturesInstallIn64BitMode=x64 diff --git a/setup.py b/setup.py index 238d996a..6610bdc5 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ setup( name="picassosr", - version="0.6.11", + version="0.7.0", author="Joerg Schnitzbauer, Maximilian T. Strauss, Rafal Kowalewski", author_email=("joschnitzbauer@gmail.com, straussmaximilian@gmail.com, rafalkowalewski998@gmail.com"), url="https://github.com/jungmannlab/picasso",