Skip to content

Commit

Permalink
WIP docs
Browse files Browse the repository at this point in the history
  • Loading branch information
lauraporta committed Oct 3, 2024
1 parent 0737e65 commit ff55b72
Show file tree
Hide file tree
Showing 3 changed files with 118 additions and 43 deletions.
22 changes: 3 additions & 19 deletions derotation/analysis/full_derotation_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -826,25 +826,9 @@ def plot_max_projection_with_center(self):
plt.savefig(self.debug_plots_folder / "max_projection_with_center.png")

def derotate_frames_line_by_line(self) -> np.ndarray:
"""Rotates the image stack line by line, using the rotation angles
by line calculated from the analog signals.
Description of the algorithm:
- takes one line from the image stack
- creates a new image with only that line
- rotates the line by the given angle
- substitutes the line in the new image
- adds the new image to the rotated image stack
Edge cases and how they are handled:
- the rotation starts in the middle of the image -> the previous lines
are copied from the first frame
- the rotation ends in the middle of the image -> the remaining lines
are copied from the last frame
Before derotation, it finds the image offset, which is the peak of
the gaussian mixture model fitted to the image histogram. It is
useful to fill in the blank pixels that appear during the derotation.
"""Wrapper for the function `derotate_an_image_array_line_by_line`.
Before calling the function, it finds the F0 image offset with
`find_image_offset`.
Returns
-------
Expand Down
117 changes: 96 additions & 21 deletions derotation/analysis/incremental_derotation_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import pandas as pd
from matplotlib.patches import Ellipse
from scipy.ndimage import rotate
from scipy.optimize import least_squares
from scipy.optimize import least_squares, OptimizeResult
from skimage.feature import blob_log
from tqdm import tqdm

Expand Down Expand Up @@ -464,23 +464,31 @@ def register_rotated_images(

return registered_images

def find_center_of_rotation(self):
def find_center_of_rotation(self) -> Tuple[int, int]:
"""Find the center of rotation by fitting an ellipse to the largest
blob centers.
Step 1: Calculate the mean images for each rotation increment.
Step 2: Find the blobs in the mean images.
Step 3: Fit an ellipse to the largest blob centers and get its center.
Returns
-------
Tuple[int, int]
The coordinates center of rotation (x, y).
"""
logging.info(
"Fitting an ellipse to the brightest blob centers "
"Fitting an ellipse to the largest blob centers "
+ "to find the center of rotation..."
)
mean_images = self.calculate_mean_images(self.image_stack)

logging.info("Finding blobs...")
(
blobs,
coord_first_blob_of_every_image,
) = self.get_coords_of_largest_blob(mean_images)
# plot blobs on top of every frame
if self.debugging_plots:
self.plot_blob_detection(blobs, mean_images)
coord_first_blob_of_every_image = self.get_coords_of_largest_blob(
mean_images
)

# Fit an ellipse to the brightest blob centers and get its center
# Fit an ellipse to the largest blob centers and get its center
center_x, center_y, a, b, theta = self.fit_ellipse_to_points(
coord_first_blob_of_every_image
)
Expand All @@ -497,8 +505,20 @@ def find_center_of_rotation(self):

return int(center_x), int(center_y)

@staticmethod
def get_coords_of_largest_blob(image_stack):
def get_coords_of_largest_blob(self, image_stack: np.ndarray) -> list:
"""Get the coordinates of the largest blob in each image.
Parameters
----------
image_stack : np.ndarray
The image stack.
Returns
-------
list
The coordinates of the largest blob in each image.
"""

blobs = [
blob_log(img, max_sigma=12, min_sigma=7, threshold=0.95, overlap=0)
for img in tqdm(image_stack)
Expand All @@ -512,17 +532,34 @@ def get_coords_of_largest_blob(image_stack):
coord_first_blob_of_every_image = [
blobs[i][0][:2].astype(int) for i in range(len(image_stack))
]

# invert x, y order
coord_first_blob_of_every_image = [
(coord[1], coord[0]) for coord in coord_first_blob_of_every_image
]

return blobs, coord_first_blob_of_every_image
# plot blobs on top of every frame
if self.debugging_plots:
self.plot_blob_detection(blobs, image_stack)

return coord_first_blob_of_every_image

def plot_blob_detection(self, blobs: list, image_stack: np.ndarray):
"""Plot the first 4 blobs in each image. This is useful to check if
the blob detection is working correctly and to see if the identity of
the largest blob is consistent across the images.
Parameters
----------
blobs : list
The list of blobs in each image.
image_stack : np.ndarray
The image stack.
"""

def plot_blob_detection(self, blobs, subgroup):
fig, ax = plt.subplots(4, 3, figsize=(10, 10))
for i, a in tqdm(enumerate(ax.flatten())):
a.imshow(subgroup[i])
a.imshow(image_stack[i])
a.set_title(f"{i*5} degrees")
a.axis("off")

Expand All @@ -535,7 +572,22 @@ def plot_blob_detection(self, blobs, subgroup):
# save the plot
plt.savefig(self.debug_plots_folder / "blobs.png")

def fit_ellipse_to_points(self, centers):
def fit_ellipse_to_points(
self, centers: list
) -> Tuple[int, int, int, int, int]:
"""Fit an ellipse to the points using least squares optimization.
Parameters
----------
centers : list
The centers of the largest blob in each image.
Returns
-------
Tuple[int, int, int, int, int]
The center of the ellipse (center_x, center_y), the semi-major
axis (a), the semi-minor axis (b), and the rotation angle (theta).
"""
# Convert centers to numpy array
centers = np.array(centers)
x = centers[:, 0]
Expand Down Expand Up @@ -580,7 +632,7 @@ def ellipse_residuals(params, x, y):
return (x_rot / a) ** 2 + (y_rot / b) ** 2 - 1

# Use least squares optimization to fit the ellipse to the points
result = least_squares(ellipse_residuals, initial_params, args=(x, y))
result: OptimizeResult = least_squares(ellipse_residuals, initial_params, args=(x, y))

# Extract optimized parameters
center_x, center_y, a, b, theta = result.x
Expand All @@ -593,8 +645,31 @@ def ellipse_residuals(params, x, y):
return center_x, center_y, a, b, theta

def plot_ellipse_fit_and_centers(
self, centers, center_x, center_y, a, b, theta
self,
centers: list,
center_x: int,
center_y: int,
a: int,
b: int,
theta: int,
):
"""Plot the fitted ellipse on the largest blob centers.
Parameters
----------
centers : list
The centers of the largest blob in each image.
center_x : int
The x-coordinate of the center of the ellipse
center_y : int
The y-coordinate of the center of the ellipse
a : int
The semi-major axis of the ellipse
b : int
The semi-minor axis of the ellipse
theta : int
The rotation angle of the ellipse
"""
# Convert centers to numpy array
centers = np.array(centers)
x = centers[:, 0]
Expand All @@ -604,7 +679,7 @@ def plot_ellipse_fit_and_centers(
fig, ax = plt.subplots(figsize=(8, 8))
# plot behind a frame of the original image
ax.imshow(self.image_stack[0], cmap="gray")
ax.scatter(x, y, label="Brightest Blob Centers", color="red")
ax.scatter(x, y, label="Largest Blob Centers", color="red")

# Plot fitted ellipse
ellipse = Ellipse(
Expand Down Expand Up @@ -636,7 +711,7 @@ def plot_ellipse_fit_and_centers(
ax.set_aspect("equal")
ax.legend()
ax.grid(True)
ax.set_title("Fitted Ellipse on brightest blob centers")
ax.set_title("Fitted Ellipse on largest blob centers")
ax.axis("off")

# save
Expand Down
22 changes: 19 additions & 3 deletions derotation/plotting_hooks/for_derotation.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,27 @@
"""
This module contains hooks for making plots while the derotation is running.
The hooks are called at specific points in the derotation process, specifically
when a line is added to the derotated image and when a frame is completed.
Also, a maximum projection plot is generated at the end of the first rotation.
They are useful for debugging purposes and for generating figures.
"""

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def line_addition(
derotated_filled_image, rotated_line, image_counter, line_counter, angle
derotated_filled_image: np.ndarray,
rotated_line: np.ndarray,
image_counter: int,
line_counter: int,
angle: float,
):
"""
Hook for plotting the derotated image and the current rotated line.
"""
fig, ax = plt.subplots(1, 2, figsize=(10, 10))

ax[0].imshow(derotated_filled_image, cmap="viridis")
Expand All @@ -26,8 +41,9 @@ def line_addition(


def image_completed(derotated_image_stack, image_counter):
"""Hook for plotting the image stack after derotation.
It is useful for debugging purposes.
"""Hook for plotting the maximum projection of the derotated image stack
after the first rotation (which ends ~149th frame) and the current
derotated image.
"""
if image_counter == 149:
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
Expand Down

0 comments on commit ff55b72

Please sign in to comment.