diff --git a/.gitignore b/.gitignore index aedc8d7..f43115e 100644 --- a/.gitignore +++ b/.gitignore @@ -81,3 +81,5 @@ venv/ # written by setuptools_scm **/_version.py + +.conda/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 24f6b0f..4e7870a 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -26,3 +26,4 @@ repos: - id: mypy additional_dependencies: - types-setuptools + - types-PyYAML diff --git a/MANIFEST.in b/MANIFEST.in index 575997a..2ae8cda 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,12 @@ include LICENSE include README.md +include *.yaml +recursive-include * *.yaml +recursive-include * *.yml +recursive-include tests *.py +recursive-include examples *.py +recursive-include tests *.png +recursive-include images *.png recursive-exclude * __pycache__ recursive-exclude * *.py[co] diff --git a/README.md b/README.md index c2fa3d8..20f02e6 100644 --- a/README.md +++ b/README.md @@ -1 +1,92 @@ -# derotation +# Derotation + +The Derotation package offers a robust solution for reconstructing rotated images, particularly those acquired during experimental preparations in the field of neuroscience research. Although it's tailored for calcium imaging data, this versatile tool is adaptable for any dataset where rotation angles are recorded, making it an essential utility for a wide array of image processing applications. + +The core algorithm, `rotate_an_image_array_line_by_line` can be also used as a standalone function to deform images by a given angle. + +## Derotation example: a grid +Derotate an image of a rotating grid using motor feedback. +On the left, the original image. On the right, the derotated image. +![Derotation example: a grid]( + ./images/rotation_by_line.png) + +## Deformation example: any image stack, any rotation angle array +With the same algorithm, you can deform any image stack by any rotation angle array. +Here an example of the deformation of a stack of images of dogs, by using a linearly increasing rotation angle array. +![Deformation example: any image stack, any rotation angle array]( + ./images/dog_rotated_by_line.png) +It has been created with the script `examples/deformation_by_line_of_a_dog.py`. + +## Introduction: Derotation in Neuroscience +In neuroscientific experiments, particularly those involving calcium imaging, precise data reconstruction is paramount. Derotation is designed to address this need, reconstructing the imaging data from experiments where the sample undergoes known rotational patterns. + +### Experimental Protocol + +The package is built to accommodate a specific experimental protocol which includes two primary recording phases: + +- **Incremental Rotation**: Samples undergo a rotation of 10 degrees every 2 seconds. This phase is crucial for sampling luminance changes due to varying positions. +- **Full Rotation**: In this key phase, the sample's activity is recorded over a complete 360-degree rotation at varying speeds and directions, following a pseudo-randomized sequence. + +## Prerequisites + +To utilize the Derotation package, the following files are required: + +- A `tif` file with the image data. +- An `aux_stim` file containing the analog signals from: + - Rotation on signal (indicating when rotation is active) + - Motor ticks (for stepper motors, one tick per angle increment) + - Frame clock (one tick per new frame acquisition) + - Line clock (one tick per new line acquisition) +- A `stumulus_randperm.mat` file detailing the stimuli randomization (including speed and direction variations). + +## Installation + +Install the Derotation package and its dependencies using pip with the following commands: + +```shell +git clone https://github.com/neuroinformatics-unit/derotation.git +cd derotation +pip install . +``` + +## Configuration +Navigate to the `derotation/config/` directory to access and edit the configuration files. You'll find two examples: `full_rotation.yml` for full rotations and `incremental_rotation.yml` for incremental rotations. + +### Setting up paths +Within these configuration files, specify the paths for your data under `paths_read` for the `tif`, `aux_stim`, and `stumulus_randperm`.mat files. The paths_write key allows you to define where the derotated TIFFs, debugging plots, and logs will be saved. + +### Config Parameters +Here's a quick rundown of the configuration parameters you'll need to adjust: + +* `channel_names`: List the names of the signals in the aux_stim file. +* `rotation_increment`: Set the motor's angle increment. +* `adjust_increment`: Enable this to adjust the rotation_increment if the motor ticks don't sum to exactly 360 degrees. +* `rot_deg`: Define a full rotation degree count. +* `debugging_plots`: Toggle this to save debugging plots. +* `analog_signals_processing`: Configure parameters for analog signal processing, including peak finding and pulse processing. +* `interpolation`: Settings related to how the rotation angle is determined from the acquisition times. + +## Usage +To run the derotation process, use one of the example scripts provided: + +```shell +python3 examples/derotate.py # For full rotation based on "full_rotation.yml" +python3 examples/derotate_incremental.py # For incremental rotation based on "incremental_rotation.yml" +``` + +### What happens behind the scenes? +**Full rotation**: +The main experimental protocol is a full rotation, where the sample is rotated 360 degrees. The rotation is recorded by the motor ticks, which are then converted to angles. The rotation angle is then interpolated to the line clock, which is used to derotate the image. Several steps are taken to ensure the accuracy of the derotation, including: +- check the number of ticks +- check the number of rotations +- adjust the rotation increment if necessary +Debugging plots can be used to visually inspect the output of each step. + +**Incremental rotation**: +The incremental rotation is a preparatory phase that precedes the full rotation. The sample is rotated 10 degrees every 2 seconds. The rotation is recorded by the motor ticks, which are then converted to angles. The rotation angle is then interpolated to the frame clock, which is used to derotate the image. The same steps are taken to ensure the accuracy of the derotation. +In addition, each frame position is corrected by using cross-correlation. + +## Explanations +### Derotation by frame vs by line +A frame-based derotation is simply the rotation of the image as a whole and is achieved using `scipy.ndimage.rotate`. This is the default derotation method for incremental rotations. +The line-based derotation is more complex and is achieved by rotating each line of the image by a certain number of pixels. This is the default derotation method for full rotations. diff --git a/derotation/analysis/full_rotation_pipeline.py b/derotation/analysis/full_rotation_pipeline.py new file mode 100644 index 0000000..0be6271 --- /dev/null +++ b/derotation/analysis/full_rotation_pipeline.py @@ -0,0 +1,882 @@ +import logging +import sys +from pathlib import Path +from typing import Tuple + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import tifffile as tiff +import yaml +from fancylog import fancylog +from scipy.signal import find_peaks +from tifffile import imsave + +from derotation.derotate_by_line import rotate_an_image_array_line_by_line +from derotation.load_data.custom_data_loaders import ( + get_analog_signals, + read_randomized_stim_table, +) + + +class FullPipeline: + """DerotationPipeline is a class that derotates an image stack + acquired with a rotating sample under a microscope. + """ + + def __init__(self, config_name): + """DerotationPipeline is a class that derotates an image stack + acquired with a rotating sample under a microscope. + In the constructor, it loads the config file, starts the logging + process, and loads the data. + It is meant to be used for the full rotation protocol, in which + the sample is rotated by 360 degrees at various speeds and + directions. + + Parameters + ---------- + config_name : str + Name of the config file without extension. + """ + self.config = self.get_config(config_name) + self.start_logging() + self.load_data() + + def __call__(self): + """Execute the steps necessary to derotate the image stack + from start to finish. + It involves: + - processing the analog signals + - rotating the image stack line by line + - adding a circular mask to the rotated image stack + - saving the masked image stack + """ + self.process_analog_signals() + rotated_images = self.rotate_frames_line_by_line() + masked = self.add_circle_mask(rotated_images) + self.save(masked) + self.save_csv_with_derotation_data() + + def get_config(self, config_name: str) -> dict: + """Loads config file from derotation/config folder. + Please edit it to change the parameters of the analysis. + + Parameters + ---------- + config_name : str + Name of the config file without extension. + Either "full_rotation" or "incremental_rotation". + + Returns + ------- + dict + Config dictionary. + """ + path_config = "derotation/config/" + config_name + ".yml" + + with open(Path(path_config), "r") as f: + config = yaml.load(f, Loader=yaml.FullLoader) + + return config + + def start_logging(self): + """Starts logging process using fancylog package. + Logs saved whenre specified in the config file. + """ + path = self.config["paths_write"]["logs_folder"] + Path(path).mkdir(parents=True, exist_ok=True) + fancylog.start_logging( + output_dir=str(path), + package=sys.modules[__name__.partition(".")[0]], + filename="derotation", + verbose=False, + ) + + def load_data(self): + """Loads data from the paths specified in the config file. + Data is stored in the class attributes. + + What is loaded: + * various parameters from config file + * image stack (tif file) + * direction and speed of rotation (from randperm file, uses \ + custom_data_loaders.read_randomized_stim_table) + * analog signals \ + (from aux file, uses `custom_data_loaders.get_analog_signals`) + + Analog signals are four files, measured in "clock_time": + * frame clock: on during acquisition of a new frame, off otherwise + * line clock: on during acquisition of a new line, off otherwise + * full rotation: when the motor is rotating + * rotation ticks: peaks at every given increment of rotation + + The data is loaded using the custom_data_loaders module, which are + sepecific to the setup used in the lab. Please edit them to load + data from your setup. + """ + logging.info("Loading data...") + + self.image_stack = tiff.imread( + self.config["paths_read"]["path_to_tif"] + ) + + self.num_frames = self.image_stack.shape[0] + self.num_lines_per_frame = self.image_stack.shape[1] + self.num_total_lines = self.num_frames * self.num_lines_per_frame + + self.direction, self.speed = read_randomized_stim_table( + self.config["paths_read"]["path_to_randperm"] + ) + + self.number_of_rotations = len(self.direction) + + ( + self.frame_clock, + self.line_clock, + self.full_rotation, + self.rotation_ticks, + ) = get_analog_signals( + self.config["paths_read"]["path_to_aux"], + self.config["channel_names"], + ) + + self.total_clock_time = len(self.frame_clock) + + self.line_use_start = self.config["interpolation"]["line_use_start"] + self.frame_use_start = self.config["interpolation"]["frame_use_start"] + + self.rotation_increment = self.config["rotation_increment"] + self.rot_deg = self.config["rot_deg"] + self.adjust_increment = self.config["adjust_increment"] + + self.filename_raw = Path( + self.config["paths_read"]["path_to_tif"] + ).stem.split(".")[0] + self.filename = self.config["paths_write"]["saving_name"] + + self.k = self.config["analog_signals_processing"]["squared_pulse_k"] + self.inter_rotation_interval_min_len = self.config[ + "analog_signals_processing" + ]["inter_rotation_interval_min_len"] + + self.debugging_plots = self.config["debugging_plots"] + + logging.info(f"Dataset {self.filename_raw} loaded") + logging.info(f"Filename: {self.filename}") + + def process_analog_signals(self): + """From the analog signals (frame clock, line clock, full rotation, + rotation ticks) calculates the rotation angles by line and frame. + + It involves: + - finding rotation ticks peaks + - identifying the rotation ticks that correspond to + clockwise and counter clockwise rotations + - removing various kinds of artifacts that derive from wrong ticks + - interpolating the angles between the ticks + - calculating the angles by line and frame + + If debugging_plots is True, it also plots: + - rotation ticks and the rotation on signal + - rotation angles by line and frame + """ + + self.rotation_ticks_peaks = self.find_rotation_peaks() + + start, end = self.get_start_end_times_with_threshold( + self.full_rotation, self.k + ) + self.rot_blocks_idx = self.correct_start_and_end_rotation_signal( + start, end + ) + self.rotation_on = self.create_signed_rotation_array() + + self.drop_ticks_outside_of_rotation() + + self.check_number_of_rotations() + if not self.is_number_of_ticks_correct() and self.adjust_increment: + ( + self.corrected_increments, + self.ticks_per_rotation, + ) = self.adjust_rotation_increment() + + self.interpolated_angles = self.get_interpolated_angles() + ( + self.line_start, + self.line_end, + ) = self.get_start_end_times_with_threshold(self.line_clock, self.k) + ( + self.frame_start, + self.frame_end, + ) = self.get_start_end_times_with_threshold(self.frame_clock, self.k) + + ( + self.rot_deg_line, + self.rot_deg_frame, + ) = self.calculate_angles_by_line_and_frame() + + if self.debugging_plots: + self.plot_rotation_on_and_ticks() + self.plot_rotation_angles() + + logging.info("✨ Analog signals processed ✨") + + def find_rotation_peaks(self) -> np.ndarray: + """Finds the peaks of the rotation ticks signal using + scipy.signal.find_peaks. It filters the peaks using + the height and distance parameters specified in the config file. + + Returns + ------- + np.ndarray + the clock times of the rotation ticks peaks + """ + + logging.info("Finding rotation ticks peaks...") + + height = self.config["analog_signals_processing"][ + "find_rotation_ticks_peaks" + ]["height"] + distance = self.config["analog_signals_processing"][ + "find_rotation_ticks_peaks" + ]["distance"] + peaks = find_peaks( + self.rotation_ticks, + height=height, + distance=distance, + )[0] + + return peaks + + @staticmethod + def get_start_end_times_with_threshold( + signal: np.ndarray, k: float + ) -> Tuple[np.ndarray, np.ndarray]: + """Finds the start and end times of the on periods of the signal. + Works for analog signals that have a squared pulse shape. + + Parameters + ---------- + signal : np.ndarray + An analog signal. + k : float + The factor used to quantify the threshold. + + Returns + ------- + Tuple(np.ndarray, np.ndarray) + The start and end times of the on periods of the signal. + """ + + mean = np.mean(signal) + std = np.std(signal) + threshold = mean + k * std + + thresholded_signal = np.zeros_like(signal) + thresholded_signal[signal > threshold] = 1 + + start = np.where(np.diff(thresholded_signal) > 0)[0] + end = np.where(np.diff(thresholded_signal) < 0)[0] + + return start, end + + def correct_start_and_end_rotation_signal( + self, + start: np.ndarray, + end: np.ndarray, + ) -> dict: + """Removes artifacts from the start and end times of the on periods + of the rotation signal. These artifacts appear as very brief off + periods that are not plausible given the experimental setup. + The two surrounding on periods are merged. + + Used the inter_rotation_interval_min_len parameter from the config + file: the minimum length of the time in between two rotations. + It is important to remove artifacts. + + Parameters + ---------- + start : np.ndarray + The start times of the on periods of rotation signal. + end : np.ndarray + The end times of the on periods of rotation signal. + + Returns + ------- + dict + Corrected start and end times of the on periods of rotation signal. + """ + + logging.info("Cleaning start and end rotation signal...") + + shifted_end = np.roll(end, 1) + mask = start - shifted_end > self.inter_rotation_interval_min_len + mask[0] = True # first rotation is always a full rotation + shifted_mask = np.roll(mask, -1) + new_start = start[mask] + new_end = end[shifted_mask] + + return {"start": new_start, "end": new_end} + + def create_signed_rotation_array(self) -> np.ndarray: + """Reconstructs an array that has the same length as the full rotation + signal. It is 0 when the motor is off, and it is 1 or -1 when the motor + is on, depending on the direction of rotation. 1 is clockwise, -1 is + counter clockwise. + Uses the start and end times of the on periods of rotation signal, and + the direction of rotation to reconstruct the array. + + Returns + ------- + np.ndarray + The rotation on signal. + """ + + logging.info("Creating signed rotation array...") + rotation_on = np.zeros(self.total_clock_time) + for i, (start, end) in enumerate( + zip( + self.rot_blocks_idx["start"], + self.rot_blocks_idx["end"], + ) + ): + rotation_on[start:end] = self.direction[i] + + return rotation_on + + def drop_ticks_outside_of_rotation(self) -> np.ndarray: + """Drops the rotation ticks that are outside of the rotation periods. + + Returns + ------- + np.ndarray + The clock times of the rotation ticks peaks, without the ticks + outside of the rotation periods. + """ + + logging.info("Dropping ticks outside of the rotation period...") + + len_before = len(self.rotation_ticks_peaks) + + rolled_starts = np.roll(self.rot_blocks_idx["start"], -1) + + # including the interval before the start of the first rotation + edited_ends = np.insert(self.rot_blocks_idx["end"], 0, 0) + rolled_starts = np.insert( + rolled_starts, 0, self.rot_blocks_idx["start"][0] + ) + + # including the end + rolled_starts[-1] = self.total_clock_time + + inter_roatation_interval = [ + idx + for i in range(self.number_of_rotations + 1) + for idx in range( + edited_ends[i], + rolled_starts[i], + ) + ] + + self.rotation_ticks_peaks = np.delete( + self.rotation_ticks_peaks, + np.where( + np.isin(self.rotation_ticks_peaks, inter_roatation_interval) + ), + ) + + len_after = len(self.rotation_ticks_peaks) + logging.info( + f"Ticks dropped: {len_before - len_after}.\n" + + f"Ticks remaining: {len_after}" + ) + + def check_number_of_rotations(self): + """Checks that the number of rotations is as expected. + + Raises + ------ + ValueError + if the number of start and end of rotations is different + ValueError + if the number of rotations is not as expected + """ + + if ( + self.rot_blocks_idx["start"].shape[0] + != self.rot_blocks_idx["end"].shape[0] + ): + raise ValueError( + "Start and end of rotations have different lengths" + ) + if self.rot_blocks_idx["start"].shape[0] != self.number_of_rotations: + raise ValueError("Number of rotations is not as expected") + + logging.info("Number of rotations is as expected") + + def is_number_of_ticks_correct(self) -> bool: + """Compares the total number of ticks with the expected number of + ticks, which is calculated from the number of rotations and the + rotation increment. + + Returns + ------- + bool + whether the number of ticks is as expected + """ + self.expected_tiks_per_rotation = ( + self.rot_deg * self.number_of_rotations / self.rotation_increment + ) + found_ticks = len(self.rotation_ticks_peaks) + + if self.expected_tiks_per_rotation == found_ticks: + logging.info(f"Number of ticks is as expected: {found_ticks}") + return True + else: + logging.warning( + f"Number of ticks is not as expected: {found_ticks}.\n" + + f"Expected ticks: {self.expected_tiks_per_rotation}\n" + + f"Delta: {found_ticks - self.expected_tiks_per_rotation}" + ) + return False + + def get_peaks_in_rotation(self, start: int, end: int) -> int: + """Counts the number of ticks in a rotation given the start and end + times of the rotation. + + Parameters + ---------- + start : int + Start clock time of the rotation. + end : int + End clock time of the rotation. + + Returns + ------- + int + The number of ticks in the rotation. + """ + return np.where( + np.logical_and( + self.rotation_ticks_peaks >= start, + self.rotation_ticks_peaks <= end, + ) + )[0].shape[0] + + def adjust_rotation_increment(self) -> Tuple[np.ndarray, np.ndarray]: + """It calculates the new rotation increment for each rotation, given + the number of ticks in each rotation. It also outputs the number of + ticks in each rotation. + + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + The new rotation increment for each rotation, and the number of + ticks in each rotation. + """ + + ticks_per_rotation = [ + self.get_peaks_in_rotation(start, end) + for start, end in zip( + self.rot_blocks_idx["start"], + self.rot_blocks_idx["end"], + ) + ] + new_increments = [self.rot_deg / t for t in ticks_per_rotation] + + logging.info(f"New increment example: {new_increments[0]:.3f}") + + return new_increments, ticks_per_rotation + + def get_interpolated_angles(self) -> np.ndarray: + """Starting from the rotation ticks and knowing the rotation increment, + it calculates the rotation angles for each clock time. + + Returns + ------- + np.ndarray + The rotation angles for each clock time. + """ + logging.info("Interpolating angles...") + + ticks_with_increment = [ + item + for i in range(self.number_of_rotations) + for item in [self.corrected_increments[i]] + * self.ticks_per_rotation[i] + ] + + cumulative_sum_to360 = np.cumsum(ticks_with_increment) % self.rot_deg + + interpolated_angles = np.zeros(self.total_clock_time) + + starts, stops = ( + self.rotation_ticks_peaks[:-1], + self.rotation_ticks_peaks[1:], + ) + + for i, (start, stop) in enumerate(zip(starts, stops)): + if cumulative_sum_to360[i + 1] < cumulative_sum_to360[i]: + cumulative_sum_to360[i] = 0 + cumulative_sum_to360[i + 1] = 0 + + interpolated_angles[start:stop] = np.linspace( + cumulative_sum_to360[i], + cumulative_sum_to360[i + 1], + stop - start, + ) + + interpolated_angles = interpolated_angles * self.rotation_on + + return interpolated_angles + + def calculate_angles_by_line_and_frame( + self, + ) -> Tuple[np.ndarray, np.ndarray]: + """From the interpolated angles, it calculates the rotation angles + by line and frame. It can use the start or the end of the line/frame to + infer the angle. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + The rotation angles by line and frame. + """ + logging.info("Calculating angles by line and frame...") + + inverted_angles = self.interpolated_angles * -1 + line_angles = np.zeros(self.num_total_lines) + frame_angles = np.zeros(self.num_frames) + + if self.line_use_start: + line_angles = inverted_angles[self.line_start] + else: + line_angles = inverted_angles[self.line_end] + + if self.frame_use_start: + frame_angles = inverted_angles[self.frame_start] + else: + frame_angles = inverted_angles[self.frame_end] + + return line_angles, frame_angles + + def clock_to_latest_line_start(self, clock_time: int) -> int: + """Get the index of the line that is being scanned at the given clock + time. + + Parameters + ---------- + clock_time : int + The clock time. + + Returns + ------- + int + The index of the line + """ + return np.where(self.line_start < clock_time)[0][-1] + + def clock_to_latest_frame_start(self, clock_time: int) -> int: + """Get the index of the frame that is being acquired at the given clock + time. + + Parameters + ---------- + clock_time : int + The clock time. + + Returns + ------- + int + The index of the frame + """ + return np.where(self.frame_start < clock_time)[0][-1] + + def clock_to_latest_rotation_start(self, clock_time: int) -> int: + """Get the index of the latest rotation that happened. + + Parameters + ---------- + clock_time : int + The clock time. + + Returns + ------- + int + The index of the latest rotation + """ + return np.where(self.rot_blocks_idx["start"] < clock_time)[0][-1] + + def plot_rotation_on_and_ticks(self): + """Plots the rotation ticks and the rotation on signal. + This plot will be saved in the debug_plots folder. + Please inspect it to check that the rotation ticks are correctly + placed during the times in which the motor is rotating. + """ + + logging.info("Plotting rotation ticks and rotation on signal...") + + fig, ax = plt.subplots(1, 1, figsize=(20, 5)) + + ax.scatter( + self.rotation_ticks_peaks, + self.rotation_on[self.rotation_ticks_peaks], + label="rotation ticks", + marker="o", + alpha=0.5, + color="orange", + ) + ax.plot( + self.rotation_on, + label="rotation on", + color="navy", + ) + + ax.plot( + self.full_rotation / np.max(self.full_rotation), + label="rotation ticks", + color="black", + alpha=0.2, + ) + + ax.set_title("Rotation ticks and rotation on signal") + + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + + plt.savefig( + Path(self.config["paths_write"]["debug_plots_folder"]) + / "rotation_ticks_and_rotation_on.png" + ) + + def plot_rotation_angles(self): + """Plots example rotation angles by line and frame for each speed. + This plot will be saved in the debug_plots folder. + Please inspect it to check that the rotation angles are correctly + calculated. + """ + logging.info("Plotting rotation angles...") + + fig, axs = plt.subplots(2, 2, figsize=(10, 10)) + + speeds = set(self.speed) + + last_idx_for_each_speed = [ + np.where(self.speed == s)[0][-1] for s in speeds + ] + last_idx_for_each_speed = sorted(last_idx_for_each_speed) + + for i, id in enumerate(last_idx_for_each_speed): + col = i // 2 + row = i % 2 + + ax = axs[col, row] + + rotation_starts = self.rot_blocks_idx["start"][id] + rotation_ends = self.rot_blocks_idx["end"][id] + + start_line_idx = self.clock_to_latest_line_start(rotation_starts) + end_line_idx = self.clock_to_latest_line_start(rotation_ends) + + start_frame_idx = self.clock_to_latest_frame_start(rotation_starts) + end_frame_idx = self.clock_to_latest_frame_start(rotation_ends) + + ax.scatter( + self.line_start[start_line_idx:end_line_idx], + self.rot_deg_line[start_line_idx:end_line_idx], + label="line angles", + color="orange", + marker="o", + ) + + ax.scatter( + self.frame_start[start_frame_idx:end_frame_idx], + self.rot_deg_frame[start_frame_idx:end_frame_idx], + label="frame angles", + color="green", + marker="o", + ) + + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + + ax.set_title( + f"Speed: {self.speed[id]}," + + f" direction: {'cw' if self.direction[id] == 1 else 'ccw'}" + ) + + fig.suptitle("Rotation angles by line and frame") + + handles, labels = ax.get_legend_handles_labels() + fig.legend(handles, labels, loc="upper right") + + plt.savefig( + Path(self.config["paths_write"]["debug_plots_folder"]) + / "rotation_angles.png" + ) + + def rotate_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 + + Returns + ------- + np.ndarray + The rotated image stack. + """ + logging.info("Starting derotation by line...") + + rotated_image_stack = rotate_an_image_array_line_by_line( + self.image_stack, + self.rot_deg_line, + ) + + logging.info("✨ Image stack rotated ✨") + return rotated_image_stack + + @staticmethod + def add_circle_mask( + image_stack: np.ndarray, + diameter: int = 256, + ) -> np.ndarray: + """Adds a circular mask to the rotated image stack. It is useful + to hide the portions of the image that are not sampled equally + during the rotation. + If a diameter is specified, the image stack is cropped to match + the diameter. The mask is then added to the cropped image stack, + and the cropped image stack is padded to match the original size. + This is important when the images are registered to correct from + motion artifacts. + + Parameters + ---------- + image_stack : np.ndarray + The image stack that you want to mask. + diameter : int, optional + The diameter of the circular mask, by default 256 + + Returns + ------- + np.ndarray + The masked image stack. + """ + img_height = image_stack.shape[1] + + # crop the image to match the new size + if diameter != img_height: + image_stack = image_stack[ + :, + int((img_height - diameter) / 2) : int( + (img_height + diameter) / 2 + ), + int((img_height - diameter) / 2) : int( + (img_height + diameter) / 2 + ), + ] + + xx, yy = np.mgrid[:diameter, :diameter] + + circle = (xx - diameter / 2) ** 2 + (yy - diameter / 2) ** 2 + mask = circle < (diameter / 2) ** 2 + + img_min = np.nanmin(image_stack) + + masked_img_array = [] + for img in image_stack: + masked_img_array.append(np.where(mask, img, img_min)) + + # pad the image to match the original size + if diameter != img_height: + delta = img_height - diameter + masked_img_array = np.pad( + masked_img_array, + ((0, 0), (int(delta / 2), int(delta / 2)), (0, 0)), + "constant", + constant_values=img_min, + ) + masked_img_array = np.pad( + masked_img_array, + ((0, 0), (0, 0), (int(delta / 2), int(delta / 2))), + "constant", + constant_values=img_min, + ) + + return np.array(masked_img_array) + + def save(self, masked: np.ndarray): + """Saves the masked image stack in the saving folder specified in the + config file. + + Parameters + ---------- + masked : np.ndarray + The masked derotated image stack. + """ + path = self.config["paths_write"]["derotated_tiff_folder"] + + imsave( + path + self.config["paths_write"]["saving_name"] + ".tif", + np.array(masked), + ) + logging.info(f"Masked image saved in {path}") + + def save_csv_with_derotation_data(self): + """Saves a csv file with the rotation angles by line and frame, + and the rotation on signal. + It is saved in the saving folder specified in the config file. + """ + df = pd.DataFrame( + columns=[ + "frame", + "rotation_angle", + "clock", + ] + ) + + df["frame"] = np.arange(self.num_frames) + df["rotation_angle"] = self.rot_deg_frame[: self.num_frames] + df["clock"] = self.frame_start[: self.num_frames] + + df["direction"] = np.nan * np.ones(len(df)) + df["speed"] = np.nan * np.ones(len(df)) + df["rotation_count"] = np.nan * np.ones(len(df)) + + rotation_counter = 0 + adding_roatation = False + for i in range(len(df)): + row = df.loc[i] + if np.abs(row["rotation_angle"]) > 0.0: + adding_roatation = True + row["direction"] = self.direction[rotation_counter] + row["speed"] = self.speed[rotation_counter] + row["rotation_count"] = rotation_counter + + df.loc[i] = row + if ( + rotation_counter < 79 + and adding_roatation + and np.abs(df.loc[i + 1, "rotation_angle"]) == 0.0 + ): + rotation_counter += 1 + adding_roatation = False + + df.to_csv( + self.config["paths_write"]["derotated_tiff_folder"] + + self.config["paths_write"]["saving_name"] + + ".csv", + index=False, + ) diff --git a/derotation/analysis/incremental_rotation_pipeline.py b/derotation/analysis/incremental_rotation_pipeline.py new file mode 100644 index 0000000..5854c6c --- /dev/null +++ b/derotation/analysis/incremental_rotation_pipeline.py @@ -0,0 +1,372 @@ +import copy +import logging +from pathlib import Path +from typing import Dict, Tuple + +import numpy as np +from matplotlib import pyplot as plt +from scipy.ndimage import rotate +from tqdm import tqdm + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +class IncrementalPipeline(FullPipeline): + """Derotate the image stack that was acquired using the incremental + rotation method. + """ + + def __init__(self, *args, **kwargs): + """Derotate the image stack that was acquired using the incremental + rotation method. + As a child of FullPipeline, it inherits all the attributes and + methods from it and adds additional logic to register the images. + Here in the constructor we specify the degrees for each incremental + rotation and the number of rotations. + """ + super().__init__(*args, **kwargs) + + self.small_rotations = 10 + self.number_of_rotations = self.rot_deg // self.small_rotations + + def __call__(self): + """Overwrite the __call__ method from the parent class to derotate + the image stack acquired using the incremental rotation method. + After processing the analog signals, the image stack is rotated by + frame and then registered using phase cross correlation. + """ + super().process_analog_signals() + rotated_images = self.roatate_by_frame() + masked_unregistered = self.add_circle_mask(rotated_images) + + mean_images = self.calculate_mean_images(masked_unregistered) + target_image = self.get_target_image(masked_unregistered) + shifts = self.get_shifts_using_phase_cross_correlation( + mean_images, target_image + ) + x_fitted, y_fitted = self.polinomial_fit(shifts) + registered_images = self.register_rotated_images( + masked_unregistered, x_fitted, y_fitted + ) + + new_diameter = self.num_lines_per_frame - 2 * max( + max(np.abs(shifts["x"])), max(np.abs(shifts["y"])) + ) + masked = self.add_circle_mask(registered_images, new_diameter) + + self.save(masked) + self.save_csv_with_derotation_data() + + def is_number_of_ticks_correct(self) -> bool: + """Overwrite the method from the parent class to check if the number + of ticks is as expected. + + Returns + ------- + bool + True if the number of ticks is as expected, False otherwise. + """ + self.expected_tiks_per_rotation = ( + self.rot_deg / self.rotation_increment + ) + found_ticks = len(self.rotation_ticks_peaks) + + if self.expected_tiks_per_rotation == found_ticks: + logging.info(f"Number of ticks is as expected: {found_ticks}") + return True + else: + logging.warning( + f"Number of ticks is not as expected: {found_ticks}.\n" + + f"Expected ticks: {self.expected_tiks_per_rotation}\n" + + f"Delta: {found_ticks - self.expected_tiks_per_rotation}" + ) + return False + + def adjust_rotation_increment(self) -> Tuple[np.ndarray, np.ndarray]: + """Overwrite the method from the parent class to adjust the rotation + increment and the number of ticks per rotation. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + The new rotation increment and the number of ticks per rotation. + """ + ticks_per_rotation = [ + self.get_peaks_in_rotation(start, end) + for start, end in zip( + self.rot_blocks_idx["start"], + self.rot_blocks_idx["end"], + ) + ] + new_increments = [self.small_rotations / t for t in ticks_per_rotation] + + logging.info(f"New increment example: {new_increments[0]:.3f}") + + return new_increments, ticks_per_rotation + + def get_interpolated_angles(self) -> np.ndarray: + """Overwrite the method from the parent class to interpolate the + angles. + + Returns + ------- + np.ndarray + The interpolated angles. + """ + logging.info("Interpolating angles...") + + ticks_with_increment = [ + item + for i in range(self.number_of_rotations) + for item in [self.corrected_increments[i]] + * self.ticks_per_rotation[i] + ] + + cumulative_sum_to360 = np.cumsum(ticks_with_increment) % self.rot_deg + + interpolated_angles = np.zeros(self.total_clock_time) + + starts, stops = ( + self.rotation_ticks_peaks[:-1], + self.rotation_ticks_peaks[1:], + ) + + for i, (start, stop) in enumerate(zip(starts, stops)): + interpolated_angles[start:stop] = np.linspace( + cumulative_sum_to360[i], + cumulative_sum_to360[i + 1], + stop - start, + ) + + return interpolated_angles * -1 + + def roatate_by_frame(self) -> np.ndarray: + """Rotate the image stack by frame. + + Returns + ------- + np.ndarray + Description of returned object. + """ + logging.info("Starting derotation by frame...") + min_value_img = np.min(self.image_stack) + new_rotated_image_stack = ( + np.ones_like(self.image_stack) * min_value_img + ) + + for idx, frame in tqdm( + enumerate(self.image_stack), total=self.num_frames + ): + rotated_img = rotate( + frame, + self.rot_deg_frame[idx], + reshape=False, + order=0, + mode="constant", + ) + rotated_img = np.where( + rotated_img == 0, min_value_img, rotated_img + ) + + new_rotated_image_stack[idx] = rotated_img + + logging.info("Finished rotating the image stack") + + return new_rotated_image_stack + + def check_number_of_frame_angles(self): + """Check if the number of rotation angles by frame is equal to the + number of frames in the image stack. + + Raises + ------ + ValueError + if the number of rotation angles by frame is not equal to the + number of frames in the image stack. + """ + if len(self.rot_deg_frame) != self.num_frames: + raise ValueError( + "Number of rotation angles by frame is not equal to the " + + "number of frames in the image stack.\n" + + f"Number of angles: {len(self.rot_deg_frame)}\n" + + f"Number of frames: {self.num_frames}" + ) + + def plot_rotation_angles(self): + """Plots example rotation angles by line and frame for each speed. + This plot will be saved in the debug_plots folder. + Please inspect it to check that the rotation angles are correctly + calculated. + """ + logging.info("Plotting rotation angles...") + + fig, ax = plt.subplots(figsize=(10, 10)) + + ax.scatter( + self.frame_start, + self.rot_deg_frame, + label="frame angles", + color="green", + marker="o", + ) + + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + + fig.suptitle("Rotation angles by frame") + + handles, labels = ax.get_legend_handles_labels() + fig.legend(handles, labels, loc="upper right") + + plt.savefig( + Path(self.config["paths_write"]["debug_plots_folder"]) + / "rotation_angles.png" + ) + + def calculate_mean_images(self, rotated_image_stack: np.ndarray) -> list: + """Calculate the mean images for each rotation angle. This required + to calculate the shifts using phase cross correlation. + + Parameters + ---------- + rotated_image_stack : np.ndarray + The rotated image stack. + + Returns + ------- + list + The list of mean images. + """ + logging.info("Calculating mean images...") + + # correct for a mismatch in the total number of frames + # and the number of angles, given by instrument error + angles_subset = copy.deepcopy(self.rot_deg_frame[2:]) + # also there is a bias on the angles + angles_subset += -0.1 + rounded_angles = np.round(angles_subset, 2) + + mean_images = [] + for i in np.arange(10, 360, 10): + images = rotated_image_stack[rounded_angles == i] + mean_image = np.mean(images, axis=0) + + mean_images.append(mean_image) + + return mean_images + + @staticmethod + def get_target_image(rotated_image_stack: np.ndarray) -> np.ndarray: + """Get the target image for phase cross correlation. This is the mean + of the first 100 images. + + Parameters + ---------- + rotated_image_stack : np.ndarray + The rotated image stack. + + Returns + ------- + np.ndarray + The target image. + """ + return np.mean(rotated_image_stack[:100], axis=0) + + def get_shifts_using_phase_cross_correlation( + self, mean_images: list, target_image: np.ndarray + ) -> Dict[str, list]: + """Get the shifts (i.e. the number of pixels that the image needs to + be shifted in order to be registered) using phase cross correlation. + + Parameters + ---------- + mean_images : list + The list of mean images for each rotation increment. + target_image : np.ndarray + The target image. + + Returns + ------- + Dict[str, list] + The shifts in x and y. + """ + + logging.info("Calculating shifts using phase cross correlation...") + shifts: Dict[str, list] = {"x": [], "y": []} + image_center = self.num_lines_per_frame / 2 + for offset_image in mean_images: + image_product = ( + np.fft.fft2(target_image) * np.fft.fft2(offset_image).conj() + ) + cc_image = np.fft.fftshift(np.fft.ifft2(image_product)) + peaks = np.unravel_index(np.argmax(cc_image), cc_image.shape) + + shift = np.asarray(peaks) - image_center + shifts["x"].append(int(shift[0])) + shifts["y"].append(int(shift[1])) + + return shifts + + def polinomial_fit(self, shifts: dict) -> Tuple[np.ndarray, np.ndarray]: + """Fit a polinomial to the shifts in order to get a smooth function + that can be used to register the images. + + Parameters + ---------- + shifts : dict + The shifts in x and y. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + The polinomial fit for x and y. + """ + logging.info("Fitting polinomial to shifts...") + + shifts["x"].insert(0, 0) + shifts["y"].insert(0, 0) + + angles_range = np.arange(0, 360, 10) + x = shifts["x"] + y = shifts["y"] + + x_fitted = np.polyfit(angles_range, x, 6) + y_fitted = np.polyfit(angles_range, y, 6) + + return x_fitted, y_fitted + + def register_rotated_images( + self, + rotated_image_stack: np.ndarray, + x_fitted: np.ndarray, + y_fitted: np.ndarray, + ) -> np.ndarray: + """Register the rotated images using the polinomial fit. + + Parameters + ---------- + rotated_image_stack : np.ndarray + The rotated image stack. + x_fitted : np.ndarray + The polinomial fit for x. + y_fitted : np.ndarray + The polinomial fit for y. + + Returns + ------- + np.ndarray + The registered image stack. + """ + + logging.info("Registering rotated images...") + registered_images = [] + for i, img in enumerate(rotated_image_stack): + angle = self.rot_deg_frame[i] + x = np.polyval(x_fitted, angle) + y = np.polyval(y_fitted, angle) + shift = (int(x), int(y)) + registered_images.append(np.roll(img, shift=shift, axis=(0, 1))) + + registered_images = np.array(registered_images) + + return registered_images diff --git a/derotation/config.yml b/derotation/config.yml deleted file mode 100644 index 7c28879..0000000 --- a/derotation/config.yml +++ /dev/null @@ -1,34 +0,0 @@ -git_version: f10323ac598209d9130983420a93506f2f6c8894 -saving: - save_to: E:\Chryssanthi\test\rot_gfap\2023_rotgfap3 - prefix: 202303271657 - suffix: 21 - index: 5 - ai_min_voltage: -10.0 - ai_max_voltage: 10.0 -nidaq: - ai: - rate: 50000.0 - channel_names: [camera,scanimage_frameclock,photodiode1,photodiode2,PI_rotCW,filtered_teensy,stage,solenoid,multiplexer_output,Vistim_ttl] - channel_ids: [ai0,ai1,ai2,ai3,ai4,ai5,ai6,ai7,ai16,ai17] - offset: [0.0,0.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.0] - scale: [1.0,1.0,1.0,1.0,1.0,40.0,-40.0,1.0,40.0,1.0] - ao: - rate: 50000.0 - channel_names: [velocity] - channel_ids: [ao0] - idle_offset: 0.5 - co: - channel_names: [camera] - channel_ids: [ctr0] - init_delay: 0 - low_samps: 107 - high_samps: 60 - clock_src: /Dev1/ai/SampleClock - do: - channel_names: [pump,multiplexer,solenoid,zero_teensy,visual_stimulus,soloist,disable_teensy] - channel_ids: [port0/line0,port0/line1,port0/line2,port0/line3,port0/line4,port0/line5,port0/line6] - clock_src: /Dev1/ai/SampleClock - di: - channel_names: [from_soloist,from_teensy] - channel_ids: [port1/line0,port1/line1] diff --git a/derotation/config/full_rotation.yml b/derotation/config/full_rotation.yml new file mode 100644 index 0000000..0ce842d --- /dev/null +++ b/derotation/config/full_rotation.yml @@ -0,0 +1,36 @@ +paths_read: + path_to_randperm: "your_path_to/stimlus_randperm.mat" + path_to_aux: "your_path_to/rotation.bin" + path_to_tif: "your_path_to/rotation.tif" +paths_write: + debug_plots_folder: "your_path_to/debug_plots/" + logs_folder: "your_path_to/logs/" + derotated_tiff_folder: "your_path_to/data_folder/" + saving_name: "derotated_image_stack" + + +channel_names: [ + "camera", + "scanimage_frameclock", + "scanimage_lineclock", + "photodiode2", + "PI_rotON", + "PI_rotticks", +] +rotation_increment: 0.2 + +adjust_increment: True +rot_deg: 360 + +debugging_plots: True + +analog_signals_processing: + find_rotation_ticks_peaks: + height: 4 + distance: 20 + squared_pulse_k: 0 + inter_rotation_interval_min_len: 1000 + +interpolation: + line_use_start: True + frame_use_start: True diff --git a/derotation/config/incremental_rotation.yml b/derotation/config/incremental_rotation.yml new file mode 100644 index 0000000..0ce842d --- /dev/null +++ b/derotation/config/incremental_rotation.yml @@ -0,0 +1,36 @@ +paths_read: + path_to_randperm: "your_path_to/stimlus_randperm.mat" + path_to_aux: "your_path_to/rotation.bin" + path_to_tif: "your_path_to/rotation.tif" +paths_write: + debug_plots_folder: "your_path_to/debug_plots/" + logs_folder: "your_path_to/logs/" + derotated_tiff_folder: "your_path_to/data_folder/" + saving_name: "derotated_image_stack" + + +channel_names: [ + "camera", + "scanimage_frameclock", + "scanimage_lineclock", + "photodiode2", + "PI_rotON", + "PI_rotticks", +] +rotation_increment: 0.2 + +adjust_increment: True +rot_deg: 360 + +debugging_plots: True + +analog_signals_processing: + find_rotation_ticks_peaks: + height: 4 + distance: 20 + squared_pulse_k: 0 + inter_rotation_interval_min_len: 1000 + +interpolation: + line_use_start: True + frame_use_start: True diff --git a/derotation/derotate.py b/derotation/derotate.py deleted file mode 100644 index 2b43536..0000000 --- a/derotation/derotate.py +++ /dev/null @@ -1,124 +0,0 @@ -from pathlib import Path - -import numpy as np -import tifffile as tiff -from matplotlib import pyplot as plt -from read_binary import read_rc2_bin -from scipy.io import loadmat -from scipy.signal import find_peaks - -# Set aux and imaging locations and initialize dip image software -rot_deg = 360 - -path = Path("/Users/lauraporta/local_data/230327_pollen") - -path_tif = path / "imaging/runtest_00001.tif" -path_aux = path / "aux_stim/202303271657_21_005.bin" -path_config = "derotation/config.yml" -path_randperm = path / "stimlus_randperm.mat" - -image = tiff.imread(path_tif) - -pseudo_random = loadmat(path_randperm) -full_rotation_blocks_direction = pseudo_random["stimulus_random"][:, 2] > 0 -dir = np.ones(4) -dir[full_rotation_blocks_direction[0:4]] = -1 - -data, dt, chan_names, config = read_rc2_bin(path_aux, path_config) -data_dict = {chan: data[:, i] for i, chan in enumerate(chan_names)} - -frame_clock = data_dict["scanimage_frameclock"] -line_clock = data_dict["camera"] -full_rotation = data_dict["PI_rotCW"] -rotation_ticks = data_dict["Vistim_ttl"] -# 0.2 deg for 1 tick - -# find the starting/ending points of the frame_clock signal -frames_start = np.where(np.diff(frame_clock) > 4)[0] -frames_end = np.where(np.diff(frame_clock) < -4)[0] -# compare if the number of frames is the same as the length of tif file -# maybe take a percentile as threshold instead of 4 - -# find the peaks of the rot_tick2 signal -rot_tick2_peaks = find_peaks(rotation_ticks, height=4)[0] -# exclude external ticks - -# identify rotation blocks -rotation_blocks = np.where(full_rotation > 4) - -# assign to each block a number in order to identify the direction - - -fig, ax = plt.subplots(4, 1, sharex=True) -ax[0].plot( - frame_clock, - label="frame clock", - color="black", - alpha=0.5, -) -# plot dots for starting and ending points of the frame_clock signal -ax[0].plot( - frames_start, - frame_clock[frames_start], - linestyle="none", - marker="o", - color="red", - alpha=0.5, -) -ax[0].plot( - frames_end, - frame_clock[frames_end], - linestyle="none", - marker="o", - color="green", - alpha=0.5, -) -ax[1].plot( - line_clock, - label="line clock", - color="red", - alpha=0.5, -) -ax[2].plot( - full_rotation, - label="rot tick", - color="blue", - alpha=0.5, -) -ax[3].plot( - rotation_ticks, - label="rot tick 2", - color="green", - alpha=0.5, - marker="o", -) -ax[3].plot( - rot_tick2_peaks, - # np.ones(len(rot_tick2_peaks)) * 5.2, - rotation_ticks[rot_tick2_peaks], - linestyle="none", - marker="*", - color="red", - alpha=0.5, -) - - -# set the initial x axis limits -for axis in ax: - # axis.set_xlim(1610000, 1800000) - # axis.set_xlim(1680000, 1710000) - axis.spines["top"].set_visible(False) - axis.spines["right"].set_visible(False) - -# set subplots titles -ax[0].set_title("Frame clock (black) and starting/ending points (red/green)") -ax[1].set_title("Line clock") -ax[2].set_title("Full rotation info") -ax[3].set_title("Rotation ticks, 0.2 deg for 1 tick (green), peaks (red)") - -# plot title -fig.suptitle("Frame clock and rotation ticks") - -plt.show() - -print("test") diff --git a/derotation/derotate_by_line.py b/derotation/derotate_by_line.py new file mode 100644 index 0000000..9500b8f --- /dev/null +++ b/derotation/derotate_by_line.py @@ -0,0 +1,110 @@ +import copy + +import numpy as np +import tqdm +from scipy.ndimage import rotate + + +def rotate_an_image_array_line_by_line( + image_stack: np.ndarray, + rot_deg_line: np.ndarray, +) -> np.ndarray: + """Rotates the image stack line by line, using the rotation angles + provided. + + 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 + + Parameters + ---------- + image_stack : np.ndarray + The image stack to be rotated. + rot_deg_line : np.ndarray + The rotation angles by line. + + Returns + ------- + np.ndarray + The rotated image stack. + """ + + num_lines_per_frame = image_stack.shape[1] + rotated_image_stack = copy.deepcopy(image_stack) + previous_image_completed = True + rotation_completed = True + + min_value_img = np.min(image_stack) + + for i, rotation in tqdm.tqdm( + enumerate(rot_deg_line), total=len(rot_deg_line) + ): + line_counter = i % num_lines_per_frame + image_counter = i // num_lines_per_frame + + is_rotating = np.absolute(rotation) > 0.00001 + image_scanning_completed = line_counter == (num_lines_per_frame - 1) + if i == 0: + rotation_just_finished = False + else: + rotation_just_finished = not is_rotating and ( + np.absolute(rot_deg_line[i - 1]) > np.absolute(rotation) + ) + + if is_rotating: + if rotation_completed and (line_counter != 0): + # when starting a new rotation in the middle of the image + rotated_filled_image = ( + np.ones_like(image_stack[image_counter]) * min_value_img + ) # non sampled pixels are set to the min val of the image + rotated_filled_image[:line_counter] = image_stack[ + image_counter + ][:line_counter] + elif previous_image_completed: + rotated_filled_image = ( + np.ones_like(image_stack[image_counter]) * min_value_img + ) + + rotation_completed = False + + img_with_new_lines = image_stack[image_counter] + line = img_with_new_lines[line_counter] + + image_with_only_line = np.zeros_like(img_with_new_lines) + image_with_only_line[line_counter] = line + + rotated_line = rotate( + image_with_only_line, + rotation, + reshape=False, + order=0, + mode="constant", + ) + + rotated_filled_image = np.where( + rotated_line == 0, rotated_filled_image, rotated_line + ) + previous_image_completed = False + if ( + image_scanning_completed and not rotation_completed + ) or rotation_just_finished: + if rotation_just_finished: + rotation_completed = True + + rotated_filled_image[line_counter + 1 :] = image_stack[ + image_counter + ][line_counter + 1 :] + + rotated_image_stack[image_counter] = rotated_filled_image + previous_image_completed = True + + return rotated_image_stack diff --git a/derotation/load_data/custom_data_loaders.py b/derotation/load_data/custom_data_loaders.py new file mode 100644 index 0000000..f22dad1 --- /dev/null +++ b/derotation/load_data/custom_data_loaders.py @@ -0,0 +1,42 @@ +import numpy as np +from scipy.io import loadmat + + +def read_randomized_stim_table(path_to_randperm): + pseudo_random = loadmat(path_to_randperm) + full_rotation_blocks_direction = pseudo_random["stimulus_random"][:, 2] > 0 + direction = np.where( + full_rotation_blocks_direction, -1, 1 + ) # 1 is counterclockwise, -1 is clockwise + + speed = pseudo_random["stimulus_random"][:, 0] + + return direction, speed + + +def read_rc2_bin(path_aux, chan_names): + n_channels = len(chan_names) + + # Read binary data saved by rc2 (int16) + # Channels along the columns, samples along the rows. + data = np.fromfile(path_aux, dtype=np.int16) + data = data.reshape((-1, n_channels)) + + # Transform the data. + # Convert 16-bit integer to volts between -10 and 10. + data = -10 + 20 * (data + 2**15) / 2**16 + + return data, chan_names + + +def get_analog_signals(path_to_aux, channel_names): + data, chan_names = read_rc2_bin(path_to_aux, channel_names) + + data_dict = {chan: data[:, i] for i, chan in enumerate(chan_names)} + + frame_clock = data_dict["scanimage_frameclock"] + line_clock = data_dict["scanimage_lineclock"] + full_rotation = data_dict["PI_rotON"] + rotation_ticks = data_dict["PI_rotticks"] + + return frame_clock, line_clock, full_rotation, rotation_ticks diff --git a/derotation/read_binary.py b/derotation/read_binary.py deleted file mode 100644 index ea37af6..0000000 --- a/derotation/read_binary.py +++ /dev/null @@ -1,31 +0,0 @@ -import numpy as np -import yaml - - -def read_rc2_bin(path_aux, path_config): - # Parse the .cfg - with open(path_config, "r") as f: - config = yaml.load(f, Loader=yaml.FullLoader) - - # Number of channels - chan_names = config["nidaq"]["ai"]["channel_names"] - n_channels = len(chan_names) - - # Read binary data saved by rc2 (int16) - # Channels along the columns, samples along the rows. - data = np.fromfile(path_aux, dtype=np.int16) - data = data.reshape((-1, n_channels)) - - # Use config file to determine offset and scale of each channel - offsets = config["nidaq"]["ai"]["offset"] - scales = config["nidaq"]["ai"]["scale"] - - # Transform the data. - # Convert 16-bit integer to volts between -10 and 10. - data = -10 + 20 * (data + 2**15) / 2**16 - # Use offset and scale to transform to correct units (cm/s etc.) - data = (data - offsets) / scales - # Also return sampling period - dt = 1 / config["nidaq"]["ai"]["rate"] - - return data, dt, chan_names, config diff --git a/docs/requirements.txt b/docs/requirements.txt index b5a754c..f6941e2 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,7 +1,10 @@ linkify-it-py myst-parser nbsphinx +numpydoc pydata-sphinx-theme setuptools-scm sphinx +sphinx-argparse sphinx-autodoc-typehints +sphinx_design diff --git a/docs/source/api_index.rst b/docs/source/api_index.rst index 348931a..e429c5a 100644 --- a/docs/source/api_index.rst +++ b/docs/source/api_index.rst @@ -1,25 +1,16 @@ API -=== +==== -math ----- +Full derotation +--------------- -.. currentmodule:: derotation.math +.. currentmodule:: derotation.analysis.full_rotation_pipeline -.. autosummary:: - :toctree: api_generated - :template: function.rst +.. autoclass:: FullPipeline - add_two_integers - subtract_two_integers +Incremental derotation +---------------------- -greetings ---------- +.. currentmodule:: derotation.analysis.incremental_rotation_pipeline -.. currentmodule:: derotation.greetings - -.. autosummary:: - :toctree: api_generated - :template: class.rst - - Greetings +.. autoclass:: IncrementalPipeline diff --git a/docs/source/conf.py b/docs/source/conf.py index 6a56fc5..ef73a02 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -25,6 +25,7 @@ author = "Laura Porta" try: release = setuptools_scm.get_version(root="../..", relative_to=__file__) + release = release.split("+")[0] except LookupError: # if git is not initialised, still allow local build # with a dummy version @@ -34,15 +35,18 @@ # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration extensions = [ - "sphinx.ext.napoleon", - "sphinx.ext.autodoc", "sphinx.ext.githubpages", - "sphinx_autodoc_typehints", + "sphinx.ext.autodoc", "sphinx.ext.autosummary", "sphinx.ext.viewcode", "sphinx.ext.intersphinx", + "sphinx.ext.napoleon", "myst_parser", + "numpydoc", "nbsphinx", + "sphinx_autodoc_typehints", + "sphinx_design", + "sphinxarg.ext", ] # Configure the myst parser to enable cool markdown features @@ -70,7 +74,24 @@ # Automatically generate stub pages for API autosummary_generate = True -autodoc_default_flags = ["members", "inherited-members"] +numpydoc_class_members_toctree = False # stops stubs warning +#toc_object_entries_show_parents = "all" +html_show_sourcelink = False + +#html_sidebars = { this is not working... +# "index": [], +# "**": [], +#} + +autodoc_default_options = { + 'members': True, + "member-order": "bysource", + 'special-members': False, + 'private-members': False, + 'inherited-members': False, + 'undoc-members': True, + 'exclude-members': "", +} # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. diff --git a/docs/source/getting_started.md b/docs/source/getting_started.md deleted file mode 100644 index bac959d..0000000 --- a/docs/source/getting_started.md +++ /dev/null @@ -1,11 +0,0 @@ -# Getting started - -Here you may demonstrate the basic functionalities your package. - -You can include code snippets using the usual Markdown syntax: - -```python -from my_package import my_function - -result = my_function() -``` diff --git a/docs/source/index.rst b/docs/source/index.rst index d180247..acfceb2 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -10,35 +10,9 @@ Welcome to derotation's documentation! :maxdepth: 2 :caption: Contents: - getting_started api_index -By default the documentation includes the following sections: -* Getting started. Here you could describe the basic functionalities of your package. To modify this page, edit the file ``docs/source/getting_started.md``. -* API: here you can find the auto-generated documentation of your package, which is based on the docstrings in your code. To modify which modules/classes/functions are included in the API documentation, edit the file ``docs/source/api_index.rst``. - -You can create additional sections with narrative documentation, -by adding new ``.md`` or ``.rst`` files to the ``docs/source`` folder. -These files shoult start with a level-1 (H1) header, -which will be used as the section title. Sub-sections can be created -with lower-level headers (H2, H3, etc.) within the same file. - -To include a section in the rendered documentation, -add it to the ``toctree`` directive in this (``docs/source/index.rst``) file. - -For example, you could create a ``docs/source/installation.md`` file -and add it to the ``toctree`` like this: - -.. code-block:: rst - - .. toctree:: - :maxdepth: 2 - :caption: Contents: - - getting_started - installation - api_index Index & Search -------------- diff --git a/examples/deformation_by_line_of_a_dog.py b/examples/deformation_by_line_of_a_dog.py new file mode 100644 index 0000000..0342f06 --- /dev/null +++ b/examples/deformation_by_line_of_a_dog.py @@ -0,0 +1,38 @@ +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +from derotation.derotate_by_line import rotate_an_image_array_line_by_line + +path = Path(__file__).parent.parent / "images/dog.png" +img = plt.imread(path) +img = img[:, :, 0] + +img_stack = np.array([[img, img, img]]).squeeze() + +img_len = img.shape[0] +rotation_angles = np.linspace(0, 180, img_len * 3) + +img_rotated = rotate_an_image_array_line_by_line(img_stack, rotation_angles) + +fig, ax = plt.subplots(1, 4, figsize=(10, 5)) + +ax[0].imshow(img, cmap="gray") +ax[0].set_title("Original image") +ax[1].imshow(img_rotated[0], cmap="gray") +ax[1].set_title("Rotated image 1") +ax[2].imshow(img_rotated[1], cmap="gray") +ax[2].set_title("Rotated image 2") +ax[3].imshow(img_rotated[2], cmap="gray") +ax[3].set_title("Rotated image 3") + +ax[0].axis("off") +ax[1].axis("off") +ax[2].axis("off") +ax[3].axis("off") + +plt.show() +plt.close() + +fig.savefig("./images/dog_rotated_by_line.png", bbox_inches="tight") diff --git a/examples/derotate.py b/examples/derotate.py new file mode 100644 index 0000000..152d4b1 --- /dev/null +++ b/examples/derotate.py @@ -0,0 +1,4 @@ +from derotation.analysis.full_rotation_pipeline import FullPipeline + +derotate = FullPipeline("full_rotation") +derotate() diff --git a/examples/derotate_incremental.py b/examples/derotate_incremental.py new file mode 100644 index 0000000..31ee7d9 --- /dev/null +++ b/examples/derotate_incremental.py @@ -0,0 +1,6 @@ +from derotation.analysis.incremental_rotation_pipeline import ( + IncrementalPipeline, +) + +derotate = IncrementalPipeline("incremental_rotation") +derotate() diff --git a/examples/derotation_slurm_job.py b/examples/derotation_slurm_job.py new file mode 100644 index 0000000..bd781ce --- /dev/null +++ b/examples/derotation_slurm_job.py @@ -0,0 +1,85 @@ +import os +import sys +from pathlib import Path + +import yaml + +from derotation.analysis.full_rotation_pipeline import FullPipeline +from derotation.analysis.incremental_rotation_pipeline import ( + IncrementalPipeline, +) + +job_id = int(sys.argv[1:][0]) +dataset_path = sys.argv[1:][1] +datasets = [path for path in os.listdir(dataset_path) if path.startswith("23")] +dataset = datasets[job_id] + +bin_files = [ + file + for file in os.listdir(f"{dataset_path}/{dataset}/aux_stim/") + if file.endswith(".bin") +] +full_rotation_bin = [file for file in bin_files if "_rotation" in file][0] +incremental_bin = [file for file in bin_files if "increment" in file][0] + +image_files = [ + file + for file in os.listdir(f"{dataset_path}/{dataset}/imaging/") + if file.endswith(".tif") +] +full_rotation_image = [file for file in image_files if "rotation_0" in file][0] +incremental_image = [file for file in image_files if "increment_0" in file][0] + +# make debug_plots and logs and derotated folder in dataset +Path(f"{dataset_path}/{dataset}/debug_plots_incremental/").mkdir( + parents=True, exist_ok=True +) +Path(f"{dataset_path}/{dataset}/debug_plots_full/").mkdir( + parents=True, exist_ok=True +) +Path(f"{dataset_path}/{dataset}/logs/").mkdir(parents=True, exist_ok=True) +Path(f"{dataset_path}/{dataset}/derotated/").mkdir(parents=True, exist_ok=True) + +for config_name in ["incremental_rotation", "full_rotation"]: + with open(f"derotation/config/{config_name}.yml") as f: + config = yaml.load(f, Loader=yaml.FullLoader) + + config["paths_read"][ + "path_to_randperm" + ] = f"{dataset_path}/stimlus_randperm.mat" + bin_name = ( + incremental_bin + if config_name == "incremental_rotation" + else full_rotation_bin + ) + config["paths_read"][ + "path_to_aux" + ] = f"{dataset_path}/{dataset}/aux_stim/{bin_name}" + image_name = ( + incremental_image + if config_name == "incremental_rotation" + else full_rotation_image + ) + config["paths_read"][ + "path_to_tif" + ] = f"{dataset_path}/{dataset}/imaging/{image_name}" + config["paths_write"][ + "debug_plots_folder" + ] = f"{dataset_path}/{dataset}/debug_plots_{config_name.split('_')[0]}" + config["paths_write"]["logs_folder"] = f"{dataset_path}/{dataset}/logs/" + config["paths_write"][ + "derotated_tiff_folder" + ] = f"{dataset_path}/{dataset}/derotated/" + config["paths_write"][ + "saving_name" + ] = f"derotated_image_stack_{config_name.split('_')[0]}" + + with open(f"derotation/config/{config_name}_{job_id}.yml", "w") as f: + yaml.dump(config, f) + + +derotate = IncrementalPipeline(f"incremental_rotation_{job_id}") +derotate() + +derotate_full = FullPipeline(f"full_rotation_{job_id}") +derotate_full() diff --git a/images/dog.png b/images/dog.png new file mode 100644 index 0000000..a9fa3d7 Binary files /dev/null and b/images/dog.png differ diff --git a/images/dog_rotated_by_line.png b/images/dog_rotated_by_line.png new file mode 100644 index 0000000..2dd6089 Binary files /dev/null and b/images/dog_rotated_by_line.png differ diff --git a/images/rotation_by_line.png b/images/rotation_by_line.png new file mode 100644 index 0000000..f629699 Binary files /dev/null and b/images/rotation_by_line.png differ diff --git a/pyproject.toml b/pyproject.toml index 6de562f..fda6607 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,18 @@ classifiers = [ "License :: OSI Approved :: BSD License", ] +dependencies = [ + "numpy", + "tiffile", + "matplotlib", + "scipy", + "PyYAML", + "fancylog", + "matplotlib", + "pandas", + "tqdm", +] + [project.urls] "Homepage" = "https://github.com/neuroinformatics-unit/derotation" "Bug Tracker" = "https://github.com/neuroinformatics-unit/derotation/issues" diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..6bd74de --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,168 @@ +import numpy as np +import pytest + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +@pytest.fixture(autouse=True) +def random(): + np.random.seed(42) + + +@pytest.fixture +def number_of_rotations(): + return 10 + + +@pytest.fixture +def rotation_len(): + return 100 + + +@pytest.fixture +def full_length(number_of_rotations, rotation_len): + return number_of_rotations * rotation_len * 3 + + +@pytest.fixture +def direction_interleaved(): + direction = np.zeros(10) + direction[::2] = 1 + direction[1::2] = -1 + return direction + + +@pytest.fixture +def direction_incremental(): + direction_1 = np.ones(5) + direction_2 = np.ones(5) * -1 + direction = np.concatenate((direction_1, direction_2)) + return direction + + +@pytest.fixture +def intervals(): + return [0, 1, 0] + + +@pytest.fixture +def dummy_signal(intervals, number_of_rotations, rotation_len): + dummy_signal = [ + num + for _ in range(number_of_rotations) + for num in intervals + for _ in range(rotation_len) + ] + return dummy_signal + + +@pytest.fixture +def full_rotation(number_of_rotations, rotation_len, intervals, dummy_signal): + dummy_signal += np.random.normal( + 0, 0.1, rotation_len * number_of_rotations * len(intervals) + ) + + return dummy_signal + + +@pytest.fixture +def start_end_times(number_of_rotations, rotation_len): + start = np.arange( + rotation_len - 1, + rotation_len - 1 + number_of_rotations * rotation_len * 3, + 3 * rotation_len, + ) + end = np.arange( + rotation_len - 1 + rotation_len, + 2 * rotation_len - 1 + number_of_rotations * rotation_len * 3, + 3 * rotation_len, + ) + + return start, end + + +@pytest.fixture +def start_end_times_with_bug(start_end_times): + start, end = start_end_times + fictional_end = 130 + fictional_start = 140 + + start = np.insert(start, 1, fictional_start) + end = np.insert(end, 0, fictional_end) + + return start, end + + +@pytest.fixture +def direction(): + direction = np.zeros(10) + direction[::2] = 1 + direction[1::2] = -1 + return direction + + +@pytest.fixture +def k(): + return 0.2 + + +@pytest.fixture +def ticks_per_rotation(): + return 100 + + +@pytest.fixture +def rotation_ticks( + full_length, + ticks_per_rotation, + number_of_rotations, +): + correct_number_of_ticks = ticks_per_rotation * number_of_rotations + number_of_ticks = correct_number_of_ticks + np.random.randint(0, 10) + + # distribute number_of_ticks in poisson indices + ticks = np.random.choice( + range(full_length), size=number_of_ticks, replace=False + ) + ticks = np.sort(ticks) + return ticks + + +@pytest.fixture +def corrected_increments(): + return np.asarray([10, 9, 11, 9, 11, 13, 14, 12, 10, 11]) + + +@pytest.fixture +def ticks_per_rotation_calculated(): + return np.asarray([36, 42, 33, 39, 34, 28, 25, 30, 36, 33]) + + +@pytest.fixture +def derotation_pipeline( + rotation_ticks, + start_end_times, + full_length, + number_of_rotations, + full_rotation, + direction, + corrected_increments, + ticks_per_rotation_calculated, + dummy_signal, +): + pipeline = FullPipeline.__new__(FullPipeline) + + pipeline.inter_rotation_interval_min_len = 50 + pipeline.rotation_ticks_peaks = rotation_ticks + pipeline.rot_blocks_idx = { + "start": start_end_times[0], + "end": start_end_times[1], + } + pipeline.number_of_rotations = number_of_rotations + pipeline.direction = direction + pipeline.total_clock_time = full_length + pipeline.full_rotation = full_rotation + pipeline.rot_deg = 360 + pipeline.rotation_on = dummy_signal + + return pipeline diff --git a/tests/test_integration/test_dropping_ticks_and_adjusting_increment.py b/tests/test_integration/test_dropping_ticks_and_adjusting_increment.py new file mode 100644 index 0000000..20d6dd4 --- /dev/null +++ b/tests/test_integration/test_dropping_ticks_and_adjusting_increment.py @@ -0,0 +1,19 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_dropping_ticks_and_adjusting_increment( + derotation_pipeline: FullPipeline, +): + len(derotation_pipeline.rotation_ticks_peaks) + derotation_pipeline.drop_ticks_outside_of_rotation() + len(derotation_pipeline.rotation_ticks_peaks) + ( + new_increments, + ticks_per_rotation, + ) = derotation_pipeline.adjust_rotation_increment() + + assert np.sum(ticks_per_rotation) == len( + derotation_pipeline.rotation_ticks_peaks + ) diff --git a/tests/test_regression/images/lenna.png b/tests/test_regression/images/lenna.png new file mode 100644 index 0000000..59ef68a Binary files /dev/null and b/tests/test_regression/images/lenna.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_1.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_1.png new file mode 100644 index 0000000..ca0b289 Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_1.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_10.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_10.png new file mode 100644 index 0000000..1d9ca8d Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_10.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_2.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_2.png new file mode 100644 index 0000000..e80f2de Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_2.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_3.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_3.png new file mode 100644 index 0000000..d61056c Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_3.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_4.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_4.png new file mode 100644 index 0000000..f9602d6 Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_4.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_5.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_5.png new file mode 100644 index 0000000..0f389da Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_5.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_6.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_6.png new file mode 100644 index 0000000..08777b1 Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_6.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_7.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_7.png new file mode 100644 index 0000000..17c8b2c Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_7.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_8.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_8.png new file mode 100644 index 0000000..48bd8e3 Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_8.png differ diff --git a/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_9.png b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_9.png new file mode 100644 index 0000000..38b8976 Binary files /dev/null and b/tests/test_regression/images/sinusoidal_rotation/rotated_lenna_9.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_1.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_1.png new file mode 100644 index 0000000..40e12e5 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_1.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_10.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_10.png new file mode 100644 index 0000000..917610f Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_10.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_2.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_2.png new file mode 100644 index 0000000..21768dd Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_2.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_3.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_3.png new file mode 100644 index 0000000..51669b8 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_3.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_4.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_4.png new file mode 100644 index 0000000..a6b10ab Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_4.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_5.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_5.png new file mode 100644 index 0000000..4e54c60 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_5.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_6.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_6.png new file mode 100644 index 0000000..53464b9 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_6.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_7.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_7.png new file mode 100644 index 0000000..1b16c49 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_7.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_8.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_8.png new file mode 100644 index 0000000..8c4d46c Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_8.png differ diff --git a/tests/test_regression/images/uniform_rotation/rotated_lenna_9.png b/tests/test_regression/images/uniform_rotation/rotated_lenna_9.png new file mode 100644 index 0000000..ef89043 Binary files /dev/null and b/tests/test_regression/images/uniform_rotation/rotated_lenna_9.png differ diff --git a/tests/test_regression/test_derotation_by_line.py b/tests/test_regression/test_derotation_by_line.py new file mode 100644 index 0000000..c7e0d28 --- /dev/null +++ b/tests/test_regression/test_derotation_by_line.py @@ -0,0 +1,65 @@ +import numpy as np +import pytest +from PIL import Image + +from derotation.analysis.full_rotation_pipeline import FullPipeline + +lenna = Image.open("tests/test_regression/images/lenna.png").convert("L") + + +@pytest.fixture +def len_stack(): + return 10 + + +@pytest.fixture +def image_stack(len_stack): + image_stack = np.array([np.array(lenna) for _ in range(len_stack)]) + return image_stack + + +@pytest.fixture +def n_lines(image_stack): + return image_stack.shape[1] + + +@pytest.fixture +def n_total_lines(image_stack, n_lines): + return image_stack.shape[0] * n_lines + + +def get_angles(kind, n_lines, n_total_lines): + if kind == "uniform": + rotation = np.linspace(0, 360, n_total_lines - n_lines) + elif kind == "sinusoidal": + rotation = ( + np.sin(np.linspace(0, 2 * np.pi, n_total_lines - n_lines)) * 360 + ) + all_angles = np.zeros(n_total_lines) + all_angles[n_lines // 2 : -n_lines // 2] = rotation + + return all_angles + + +def test_rotation_by_line(image_stack, n_lines, n_total_lines, len_stack): + pipeline = FullPipeline.__new__(FullPipeline) + pipeline.image_stack = image_stack + + for kind in ["uniform", "sinusoidal"]: + pipeline.rot_deg_line = get_angles(kind, n_lines, n_total_lines) + pipeline.num_lines_per_frame = n_lines + + rotated_images = pipeline.rotate_frames_line_by_line() + + assert len(rotated_images) == len_stack + assert rotated_images[0].shape == (n_lines, n_lines) + + for i, image in enumerate(rotated_images): + target_image = Image.open( + "tests/test_regression/images/" + + f"{kind}_rotation/rotated_lenna_{i + 1}.png" + ) + target_image = np.array(target_image.convert("L")) + assert np.allclose( + image, target_image, atol=1 + ), f"Failed for {kind} rotation, image {i + 1}" diff --git a/tests/test_unit/test_adjust_rotation_increment.py b/tests/test_unit/test_adjust_rotation_increment.py new file mode 100644 index 0000000..318617a --- /dev/null +++ b/tests/test_unit/test_adjust_rotation_increment.py @@ -0,0 +1,42 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_adjust_rotation_increment_360( + derotation_pipeline: FullPipeline, + corrected_increments, + ticks_per_rotation_calculated, +): + ( + new_increments, + new_ticks_per_rotation, + ) = derotation_pipeline.adjust_rotation_increment() + + new_increments = np.round(new_increments, 0) + + assert np.all(new_increments == corrected_increments), f"{new_increments}" + assert np.all( + new_ticks_per_rotation == ticks_per_rotation_calculated + ), f"{new_ticks_per_rotation}" + + +def test_adjust_rotation_increment_5( + derotation_pipeline: FullPipeline, +): + derotation_pipeline.rot_deg = 5 + + ( + new_increments, + _, + ) = derotation_pipeline.adjust_rotation_increment() + + new_increments = np.round(new_increments, 3) + + correct_increments = np.array( + [0.139, 0.119, 0.152, 0.128, 0.147, 0.179, 0.2, 0.167, 0.139, 0.152] + ) + + assert np.all( + new_increments == correct_increments + ), f"new_increments: {new_increments}" diff --git a/tests/test_unit/test_create_signed_rotation_array.py b/tests/test_unit/test_create_signed_rotation_array.py new file mode 100644 index 0000000..bfe8ae7 --- /dev/null +++ b/tests/test_unit/test_create_signed_rotation_array.py @@ -0,0 +1,29 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_create_signed_rotation_array_interleaved( + derotation_pipeline: FullPipeline, + start_end_times: tuple, +): + start, end = start_end_times + rotation_on = derotation_pipeline.create_signed_rotation_array() + + for idx in range(0, len(start), 2): + assert np.all(rotation_on[start[idx] : end[idx]] == 1) + assert np.all(rotation_on[start[idx + 1] : end[idx + 1]] == -1) + + +def test_create_signed_rotation_array_incremental( + derotation_pipeline: FullPipeline, + start_end_times: tuple, + direction_incremental: np.ndarray, +): + derotation_pipeline.direction = direction_incremental + start, end = start_end_times + rotation_on = derotation_pipeline.create_signed_rotation_array() + + for idx in range(0, 5): + assert np.all(rotation_on[start[idx] : end[idx]] == 1) + assert np.all(rotation_on[start[idx + 5] : end[idx + 5]] == -1) diff --git a/tests/test_unit/test_drop_ticks.py b/tests/test_unit/test_drop_ticks.py new file mode 100644 index 0000000..c55e8c6 --- /dev/null +++ b/tests/test_unit/test_drop_ticks.py @@ -0,0 +1,12 @@ +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_drop_ticks_generated_randomly( + derotation_pipeline: FullPipeline, +): + len_before = len(derotation_pipeline.rotation_ticks_peaks) + derotation_pipeline.drop_ticks_outside_of_rotation() + len_after = len(derotation_pipeline.rotation_ticks_peaks) + + assert len_before > len_after + assert len(derotation_pipeline.rotation_ticks_peaks) == 333 diff --git a/tests/test_unit/test_finding_correct_start_end_times_with_threshold.py b/tests/test_unit/test_finding_correct_start_end_times_with_threshold.py new file mode 100644 index 0000000..7093245 --- /dev/null +++ b/tests/test_unit/test_finding_correct_start_end_times_with_threshold.py @@ -0,0 +1,21 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_finding_correct_start_end_times_with_threshold( + derotation_pipeline: FullPipeline, + full_rotation: np.ndarray, + k: int, + number_of_rotations: int, + rotation_len: int, +): + start, end = derotation_pipeline.get_start_end_times_with_threshold( + full_rotation, k + ) + + assert len(start) == len(end) + assert start[0] < end[0] + assert end[-1] > start[-1] + assert len(start) == number_of_rotations + assert start[0] == rotation_len - 1 diff --git a/tests/test_unit/test_get_interpolated_angles.py b/tests/test_unit/test_get_interpolated_angles.py new file mode 100644 index 0000000..ebcd070 --- /dev/null +++ b/tests/test_unit/test_get_interpolated_angles.py @@ -0,0 +1,30 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_get_interpolated_angles(derotation_pipeline: FullPipeline): + derotation_pipeline.drop_ticks_outside_of_rotation() + ( + derotation_pipeline.corrected_increments, + derotation_pipeline.ticks_per_rotation, + ) = derotation_pipeline.adjust_rotation_increment() + + angles = derotation_pipeline.get_interpolated_angles() + angles = np.round(angles, 0) + + test_angles = [ + 0, + 10, + 21, + 21, + 31, + 31, + 34, + 38, + 41, + 41, + 51, + ] + + assert np.all(angles[99:110] == test_angles) diff --git a/tests/test_unit/test_placeholder.py b/tests/test_unit/test_placeholder.py deleted file mode 100644 index 3ada1ee..0000000 --- a/tests/test_unit/test_placeholder.py +++ /dev/null @@ -1,2 +0,0 @@ -def test_placeholder(): - assert True diff --git a/tests/test_unit/test_removing_brief_off_periods.py b/tests/test_unit/test_removing_brief_off_periods.py new file mode 100644 index 0000000..4f7e961 --- /dev/null +++ b/tests/test_unit/test_removing_brief_off_periods.py @@ -0,0 +1,27 @@ +import numpy as np + +from derotation.analysis.full_rotation_pipeline import FullPipeline + + +def test_removing_brief_off_periods( + start_end_times: tuple, + start_end_times_with_bug: tuple, + derotation_pipeline: FullPipeline, +): + start_buggy, end_buggy = start_end_times_with_bug + start, end = start_end_times + + corrected = derotation_pipeline.correct_start_and_end_rotation_signal( + start_buggy, end_buggy + ) + + assert len(corrected["start"]) == len(corrected["end"]) + assert len(corrected["start"]) == len(start) + assert len(corrected["end"]) == len(end) + assert corrected["start"][0] == start[0] + assert corrected["end"][0] == end[0] + + assert np.any( + corrected["end"] - corrected["start"] + >= derotation_pipeline.inter_rotation_interval_min_len + )