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/.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} diff --git a/changelog.rst b/changelog.rst index 16e9f796..9eb18fe7 100644 --- a/changelog.rst +++ b/changelog.rst @@ -1,7 +1,16 @@ Changelog ========= -Last change: 06-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/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/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/__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/picasso/aim.py b/picasso/aim.py new file mode 100644 index 00000000..5939e591 --- /dev/null +++ b/picasso/aim.py @@ -0,0 +1,728 @@ +""" + 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 + + +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, seg_bounds, intersect_d, roi_r, width, + aim_round=1, 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. + 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 + 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. + 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. + + 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. + """ + + assert aim_round in [1, 2], "aim_round must be 1 or 2." + + # number of segments + n_segments = len(seg_bounds) - 1 + 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 + start_idx = 1 if aim_round == 1 else 0 + if progress is not None: + iterator = range(start_idx, n_segments) + else: + 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 > 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] + + # 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 + + # 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 + 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(seg_bounds[-1]) + 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, 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 + 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 = len(seg_bounds) - 1 + 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 + start_idx = 1 if aim_round == 1 else 0 + if progress is not None: + iterator = range(start_idx, n_segments) + else: + 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 > 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] + + # 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 + + # 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 + t = (seg_bounds[1:] + seg_bounds[:-1]) / 2 + drift_z_pol = _InterpolatedUnivariateSpline(t, drift_z, k=3) + t_inter = _np.arange(seg_bounds[-1]) + 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"] + n_frames = info[0]["Frames"] + + # frames should start at 1 + frame = locs["frame"] + 1 + # 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] + 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, 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, 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 + 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, 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, 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) + 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..db8af49f 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 segment (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"] @@ -9682,7 +9795,6 @@ def undrift(self): "Correlating image pairs", 0, n_pairs, self ) try: - start_time = time.time() # find drift and apply it to locs drift, _ = postprocess.undrift( locs, @@ -9692,11 +9804,6 @@ def undrift(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 @@ -9721,23 +9828,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 +9957,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 +10981,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 +10997,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/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) 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/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 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/picasso/localize.py b/picasso/localize.py index 89de078c..bc3c0b23 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) @@ -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): 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 1374679a..ea96bea2 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. Previous versions ----------------- @@ -126,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 -------------- @@ -136,6 +132,13 @@ 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 `__ +- 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 ------- 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",