1. Place a power meter sensor (e.g., Thorlabs S170C) on the microscope stage.

2. Center the sensor with the objective to ensure proper alignment.

3. Zero the sensor to ensure accurate readings.

4. Using your power meter controller (e.g., Thorlabs PM400) or compatible software, select the wavelength of the light source you want to monitor.

5. Turn on the light source and immediately record the power output over time until it stabilizes.

   I personally record every 10 seconds for 24 hours and stops when it has been stable for 1 h.

6. Repeat steps 3 to 5 for each light source you want to monitor

   Keep the light source turned on at all times. Depending on your hardware, the light source may remain continuously on or be automatically shut down by the software when not in use.

Use the provided spreadsheet template, Illumination_Warmup Kinetic_Template.xlsx, and fill in the orange cells to visualize your results. For each light source, plot the measured power output (in mW) against time to analyze the data.

Calculate the relative power using the formula: Relative Power = (Power / Max Power). Then, plot the Relative Power (%) against time to visualize the data.

We observe some variability in the power output for the 385 nm light source.

To assess the stability:

1. Define a Stability Duration Window:
   Select a time period (e.g., 10 minutes) during which the power output should remain stable.

2. Specify a Maximum Coefficient of Variation (CV) Threshold:
   Determine an acceptable variability limit for the selected window (e.g., 0.01%).

3. Calculate the Coefficient of Variation (CV):
   Use the formula: CV = (Standard Deviation / Mean)
   Compute the CV for the specified stability duration window.

4. Visualize Stability:
   Plot the calculated CV over time to analyze the stability of the power output.

We observe that the light sources stabilize quickly, within less than 10 minutes, while the 385 nm light source takes approximately 41 minutes to reach stability.

Report the results in a table

 Selected Stability Duration Window (min): 10 min and Maximum Coefficient of Variation: 0.01%

The illumination warm-up time for this instrument is approximately 40 minutes. This duration is essential for ensuring accurate quantitative measurements, as the Coefficient of Variation (CV) threshold is strict, with a maximum allowable variation of 0.01% within a 10-minute window.

1. Warm up the light sources for the required duration (see previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Center the sensor with the objective to ensure proper alignment.
4. Zero the sensor to ensure accurate readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn the light source on to 100% power.
7. Record the average power output over 10 seconds.
   I personally re-use the data collected during the warm-up kinetic experiment (when the power is stable) for this purpose.
8. Repeat steps 5 to 7 for each light source and wavelength you wish to monitor.

Fill in the orange cells in the Illumination_Maximum Power Output_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured maximum power output (in mW).

Plot the measured maximum power output (in mW) and compare it to the manufacturer’s specifications. Calculate the Relative Power using the formula: Relative Power = (Measured Power / Manufacturer Specifications) and plot the relative power for each light source

Report the results in a table

This instrument provides 80% of the power specified by the manufacturer. These results are consistent, as the manufacturer’s specifications are based on a different objective, and likely different filters and mirrors, which can affect the measured power output.

1. Warm up the light sources for the required duration (see previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Center the sensor with the objective to ensure proper alignment.
4. Zero the sensor to ensure accurate readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn the light source on to 100% power.
7. Record the power output every 100 ms for 1 minute. For some microscope dedicated to a fast imaging it might be required to record stability at a faster rate. The THorlab S170C sensor can record at 1 kHz
   

   I personally acquire this data immediately after the warm-up kinetic experiment, without turning off the light source.

8. Repeat steps 5 to 7 for each light source and wavelength you wish to monitor.

Fill in the orange cells in the Illumination_Stability_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured power output (in mW) over time.

Calculate the relative power using the formula: Relative Power = (Power / Max Power). Then, plot the Relative Power (%) over time.

Calculate the Stability Factor (S) using the formula: S (%) = 100 × (1 - (Pmax - Pmin) / (Pmax + Pmin)). Also, calculate the Coefficient of Variation (CV) using the formula: CV = Standard Deviation / Mean. Reports the results in a table.

From the Stability Factor results, we observe that the difference between the maximum and minimum power is less than 0.03%. Additionally, the Coefficient of Variation indicates that the standard deviation is less than 0.004% of the mean value, demonstrating excellent power stability.

The light sources exhibit a very high stability, >99.9% during a 1-minute period.

1. Warm up the light sources for the required duration (see previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Ensure the sensor is properly centered with the objective for accurate measurements.
4. Zero the sensor to ensure precise readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the desired wavelength of the light source.
6. Turn on the light source to 100% intensity.
7. Record the power output every 10 seconds for 15 minutes.

   I personally re-use the data collected during the warm-up kinetic experiment (when the power is stable) for this purpose.

8. Repeat steps 5 to 7 for each light source you wish to monitor.

Fill in the orange cells in the Illumination_Stability_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured power output (in mW) over time.

Calculate the relative power using the formula: Relative Power = Power / Max Power. Then, plot the Relative Power (%) over time.

Calculate the Stability Factor (S) using the formula: S (%) = 100 × (1 - (Pmax - Pmin) / (Pmax + Pmin)). Also, calculate the Coefficient of Variation (CV) using the formula: CV = Standard Deviation / Mean. Reports the results in a table.

From the Stability Factor results, we observe that the difference between the maximum and minimum power is less than 0.01%. Additionally, the Coefficient of Variation indicates that the standard deviation is less than 0.004% of the mean value, demonstrating excellent power stability.

The light sources exhibit high stability, maintaining >99.9% stability during a 15-minute period.

1. Warm up the light sources for the required duration (see previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Ensure the sensor is properly centered with the objective for accurate measurements.
4. Zero the sensor to ensure precise readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the desired wavelength of the light source.
6. Turn on the light source to 100% intensity.
7. Record the power output every 10 seconds for 1 hour.

   I personally re-use the data collected during the warmup kinetic experiment.

8. Repeat steps 5 to 7 for each light source you wish to monitor

Fill in the orange cells in the Illumination_Stability_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured power output (in mW) over time to assess stability.

Calculate the relative power using the formula: Relative Power = Power / Max Power. Then, plot the Relative Power (%) over time.

Calculate the Stability Factor (S) using the formula: S (%) = 100 × (1 - (Pmax - Pmin) / (Pmax + Pmin)). Also, calculate the Coefficient of Variation (CV) using the formula: CV = Standard Deviation / Mean. Reports the results in a table.

From the Stability Factor results, we observe that the difference between the maximum and minimum power is less than 0.03%. Additionally, the Coefficient of Variation indicates that the standard deviation is less than 0.02% of the mean value, demonstrating excellent power stability.

The light sources exhibit exceptional stability, maintaining a performance of >99.9% during a 1-hour period.

1. Warm up the light sources for the required duration (see previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Center the sensor with the objective to ensure proper alignment.
4. Zero the sensor to ensure accurate readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn the light source on to 100% power.
7. Record the average power output over 10 seconds.
   I personally re-use the data collected for the maximal power output section and plot it over time.
8. Repeat steps 5 to 7 for each light source and wavelength you wish to monitor.

Fill in the orange cells in the Illumination_Stability_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured power output (in mW) over time to assess the stability of the illumination.

Calculate the relative power using the formula: Relative Power = Power / Max Power. Then, plot the Relative Power (%) over time.

Calculate the Relative Power_Spec by comparing the measured power to the manufacturer’s specifications using the following formula: Relative Power_Spec = Power / Power_Spec. Then, plot the Relative Power_Spec (% Spec) over time.

We expect a gradual decrease in power output over time due to the aging of hardware. Light sources should be replaced when the Relative Power_Spec falls below 50%.

Reports the results in a table.

The light sources are somehow stable over the last 2 years but a decrease in the maximum power output is seen.

1. Warm up the light sources for the required duration (refer to the previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Ensure the sensor is properly centered with the objective for accurate measurements.
4. Zero the sensor to ensure precise readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn on the light source and adjust its intensity to 0%, then incrementally increase to 10%, 20%, 30%, and so on, up to 100%.
   I typically collect this data immediately after the warm-up kinetic phase and once the real-time power stability data has been recorded.
7. Record the power output corresponding to each input level.
8. Repeat steps 5 to 7 for each light source you wish to monitor

Fill in the orange cells in the Illumination_Linearity_Template.xlsx spreadsheet to visualize your results. For each light source, plot the measured power output (in mW) as a function of the input (%).

Calculate the Relative Power using the formula: Relative Power = Power / MaxPower. Then, plot the Relative Power (%) as a function of the input (%).

Determine the equation for each curve, which is typically a linear relationship of the form: Output = K × Input. Report the slope (K) and the coefficient of determination (R²) for each curve in a table.

The slopes demonstrate a nearly perfect linear relationship between the input and the measured output power, with values very close to 1. The coefficient of determination (R²) indicates a perfect linear fit, showing no deviation from the expected relationship.

The light sources are highly linear.

1. Warm up the light sources for the required duration (refer to the previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Ensure the sensor is properly centered with the objective.
4. Zero the sensor to ensure accurate readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn on the light source to 100% intensity.
7. Record the power output for each objective, as well as without the objective in place.

   I typically collect this data after completing the warm-up kinetic phase, followed by the real-time power stability measurements, and immediately after recording the power output linearity.

8. Repeat steps 5 to 7 for each light source and wavelength

Fill in the orange cells in the Objective and cube transmittance_Template.xlsx spreadsheet to visualize your results. For each objective, plot the measured power output (in mW) as a function of the wavelength (in nm).

Calculate the Relative Transmittance using the formula: Relative Transmittance = Power / PowerNoObjective. Then, plot the Relative Transmittance (%) as a function of the wavelength (in nm).

Calculate the average transmittance for each objective and report the results in a table. Compare the average transmittance to the specifications provided by the manufacturer to assess performance.

The measurements are generally close to the specifications, with the exception of the 63x-1.4 objective. This deviation is expected, as the 63x objective has a smaller back aperture, which reduces the amount of light it can receive. Additionally, you can compare the shape of the transmittance curves to further assess performance.

The objectives are transmitting light properly.

1. Warm up the light sources for the required duration (refer to the previous section).
2. Place the power meter sensor (e.g., Thorlabs S170C) on the microscope stage.
3. Ensure the sensor is properly centered with the objective.
4. Zero the sensor to ensure accurate readings.
5. Using your power meter controller (e.g., Thorlabs PM400) or software, select the wavelength of the light source you wish to monitor.
6. Turn on the light source to 100% intensity.
7. Record the power output for each filter cube.

   I typically collect this data after completing the warm-up kinetic phase, followed by the real-time power stability measurements, and immediately after recording the power output linearity.

8. Repeat steps 5 to 7 for each light source and wavelength

Fill in the orange cells in the Objective and cube transmittance_Template.xlsx spreadsheet to visualize your results. For each filter cube, plot the measured power output (in mW) as a function of the wavelength (in nm).

Calculate the Relative Transmittance using the formula: Relative Transmittance = Power / PowerMaxFilter. Then, plot the Relative Transmittance (%) as a function of the wavelength (in nm).

Calculate the average transmittance for each filter at the appropriate wavelengths and report the results in a table.

• The DAPI cube transmits only 14% of the excitation light compared to the Quad Band Pass DAPI/GFP/Cy3/Cy5. While it is still usable, it will provide a low signal. This is likely because the excitation filter within the cube does not match the light source properly. Since an excitation filter is already included in the light source, the filter in this cube could be removed.

• The GFP and DsRed cubes transmit 47% of the excitation light compared to the Quad Band Pass DAPI/GFP/Cy3/Cy5, and they are functioning properly.

• The DHE cube does not transmit any light from the Colibri. This cube may need to be removed and stored.

• The Cy5 cube transmits 84% of the excitation light compared to the Quad Band Pass DAPI/GFP/Cy3/Cy5, and it is working properly.

Actions to be Taken:

• Remove the excitation filter from the DAPI cube, as it does not properly match the light source and is redundant with the excitation filter already included in the light source.
• Remove and store the DHE cube, as it does not transmit any light from the Colibri and is no longer functional.

1. Place a fluorescent plastic slide (Thorlabs FSK5) onto the stage.
2. Center the slide and the objective.
3. Adjust the focus to align with the surface of the slide.
   

   I typically use a red fluorescent slide and focus on a scratch mark on its surface for alignment.

4. Slightly adjust the focus deeper into the slide to minimize the visibility of dust, dirt, and scratches.

5. Modify the acquisition parameters to ensure the image is properly exposed.

   I typically use the auto-exposure feature, aiming for a targeted intensity of 30%.

6. Capture a multi-channel image.
7. Repeat steps 5 and 6 for each objective and filter combination.

You should have acquired several multi-channel images that now need processing to yield meaningful results. To process them, use the Field Illumination analysis feature of the MetroloJ_QC plugin for FIJI. For more information about the MetroloJ_QC plugin please refer to manual available at the MontpellierRessourcesImagerie on GitHub.

1. Open FIJI.
2. Load your image by dragging it into the FIJI bar.
3. Launch the MetroloJ_QC plugin by navigating to Plugins > MetroloJ QC.
4. Click on Field Illumination Report.
5. Enter a title for your report.
6. Type in your name.
7. Click Filter Parameters and input the filter's names, excitation, and emission wavelengths.

8. Check Remove Noise using Gaussian Blur.

9. Enable Apply Tolerance to the Report and reject uniformity and accuracy values below 80%.

10. Click File Save Options.

11. Select Save Results as PDF Reports

12. Select Save Results as spreadsheets.

13. Click OK.

14. Repeat steps 4 through 13 for each image you have acquired.

This will generate detailed results stored in a folder named Processed. The Processed folder will be located in the same directory as the original images, with each report and result saved in its own sub-folder.

• Use the following R script Process Field Uniformity Results.R to merge all _results.xls files into a single file.

• This script will load all _results.xls files located in the "Input" folder on your desktop or prompt you to select an input folder.

• It will compile and merge data generated by the MetroloJ_QC Field Illumination analysis, including metrics such as Uniformity and Centering Accuracy for each channel from every file.

• Additionally, it will split the filenames using the _ character and create columns named Variable-001, Variable-002, etc. This will help in organizing the data if the filenames are formatted as Objective_VariableA_VariableB, for example.

• The combined data will be saved as Combined_output.csv in the Output folder on your desktop.

# Clear the workspace
rm(list = ls())

# Set default input and output directories
default_input_dir <- file.path(Sys.getenv("USERPROFILE"), "Desktop", "Input")  # Default to "Input" on Desktop
InputFolder <- default_input_dir  # Use default folder

# Prompt user to select a folder (optional, currently commented)
InputFolder <- choose.dir(default = default_input_dir, caption = "Select an Input Folder") # You may comment this line if you want to directly use the InputFolder on your desktop

# If no folder is selected, fall back to the default
if (is.na(InputFolder)) {
  InputFolder <- default_input_dir
}

# Specify the Output folder (this is a fixed folder on the Desktop)
OutputFolder <- file.path(Sys.getenv("USERPROFILE"), "Desktop", "Output")
if (!dir.exists(OutputFolder)) dir.create(OutputFolder, recursive = TRUE)

# List all XLS files ending with "_results.xls" in the folder and its subfolders
xls_files <- list.files(path = InputFolder, pattern = "_results.xls$", full.names = TRUE, recursive = TRUE)
if (length(xls_files) == 0) stop("No _results.xls files found in the selected directory.")

# Define the header pattern to search for
header_pattern <- "Filters combination/set.*Uniformity.*Field Uniformity.*Centering Accuracy.*Image.*Coef. of Variation.*Mean c fit value"

# Initialize a list to store tables
all_tables <- list()

# Loop over all XLS files
for (file in xls_files) {
  # Read the file as raw lines (since it's tab-separated)
  lines <- readLines(file)
  
  # Find the line that matches the header pattern
  header_line_index <- grep(header_pattern, lines)
  
  if (length(header_line_index) == 0) {
    cat("No matching header found in file:", file, "\n")
    next  # Skip this file if the header isn't found
  }
  
  # Extract the table starting from the header row
  start_row <- header_line_index[1]
  
  # The table might end at an empty line, so let's find the next empty line after the header
  end_row <- min(grep("^\\s*$", lines[start_row:length(lines)]), na.rm = TRUE) + start_row - 1
  if (is.na(end_row)) end_row <- length(lines)  # If no empty line found, use the last line
  
  # Extract the table rows (lines between start_row and end_row)
  table_lines <- lines[start_row:end_row]
  
  # Convert the table lines into a data frame
  table_data <- read.table(text = table_lines, sep = "\t", header = TRUE, stringsAsFactors = FALSE)
  
  # Add the filename as the first column (without extension)
  filename_no_extension <- tools::file_path_sans_ext(basename(file))
  table_data$Filename <- filename_no_extension
  
  # Split the filename by underscores and add as new columns
  filename_parts <- unlist(strsplit(filename_no_extension, "_"))
  for (i in 1:length(filename_parts)) {
    table_data[[paste0("Variable-", sprintf("%03d", i))]] <- filename_parts[i]
  }
  
  # Rename columns based on the exact names found in the data
  colnames(table_data) <- gsub("Filters.combination.set", "Channel", colnames(table_data))
  colnames(table_data) <- gsub("Uniformity....", "Uniformity (%)", colnames(table_data))
  colnames(table_data) <- gsub("Field.Uniformity....", "Field Uniformity (%)", colnames(table_data))
  colnames(table_data) <- gsub("Centering.Accuracy....", "Centering Accuracy (%)", colnames(table_data))
  colnames(table_data) <- gsub("Coef..of.Variation", "CV", colnames(table_data))
  
  # Remove the "Mean.c.fit.value" column if present
  colnames(table_data) <- gsub("Mean.c.fit.value", "", colnames(table_data))
  
  # Drop any columns that were renamed to an empty string
  table_data <- table_data[, !grepl("^$", colnames(table_data))]
  
  # Divide the relevant columns by 100 to convert percentages into decimal values
  table_data$`Uniformity (%)` <- table_data$`Uniformity (%)` / 100
  table_data$`Field Uniformity (%)` <- table_data$`Field Uniformity (%)` / 100
  table_data$`Centering Accuracy (%)` <- table_data$`Centering Accuracy (%)` / 100
  
  # Store the table in the list with the file name as the name
  all_tables[[basename(file)]] <- table_data
}

# Check if any tables were extracted
if (length(all_tables) == 0) {
  stop("No tables were extracted from the files.")
}

# Optionally, combine all tables into a single data frame (if the structures are consistent)
combined_data <- do.call(rbind, all_tables)

# View the combined data
#head(combined_data)

# Optionally, save the combined data to a CSV file
output_file <- file.path(OutputFolder, "Field-Illumination_Combined-Data.csv")
write.csv(combined_data, output_file, row.names = FALSE)
cat("Combined data saved to:", output_file, "\n")

The following spreadsheet provides a dataset that can be manipulated with a pivot table to generate informative graphs and statistics Illumination_Uniformity_Template.xlsx. Plot the uniformity and centering accuracy for each objective.

 Plot the uniformity and centering accuracy for each filter set. 

This specific instrument has a quad-band filter as well as individual filter cubes. We can plot the uniformity and centering accuracy per filter types. 

The uniformity and centering accuracy are excellent across all objectives and filters, consistently exceeding 90%. However, the single-band filter cubes exhibit slightly lower uniformity and centering accuracy compared to the quad-band filter cube.

1. Place 4 µm diameter fluorescent beads (TetraSpec Fluorescent Microspheres Size Kit, mounted on a slide) on the microscope stage.

2. Center an isolated bead under a high-NA dry objective.

3. Crop the acquisition area to only visualize one bead but keep it large enough to anticipate a potential drift along the X and Y axes (100 um FOV should be enough)
   

4. Select an imaging channel appropriate for the fluorescent beads

   I typically use the Cy5 channel which is very bright and resistant to bleaching, yet this channel has a lower resolution, but it does not really matter here

5. Acquire a large Z-stack at 1-minute intervals for a duration of 24 hours.

   To ensure accurate measurements, it is essential to account for potential drift along the Z-axis by acquiring a Z-stack that is substantially larger than the visible bead size. I typically acquire a 40 µm Z-stack.
   

1. Open your image in FIJI
2. If necessary, crop the image to focus on a single bead for better visualization.
3. Use the TrackMate plugin included in FIJI to detect and track spots over time Plugins\Tracking\TrackMate.
4. Apply Difference of Gaussians (DoG) spot detection with a detection size of 4 µm
5. Enable sub-pixel localization for increased accuracy
6. Click Preview to visualize the spot detection
7. Set a quality threshold (click and slide) high enough to detect a single spot per frame
8. Click Next and follow the detection and tracking process
9. Save the detected Spots coordinates as a CSV file for further analysis

Open the spreadsheet template XYZ Drift Kinetic_Template.xlsx and fill in the orange cells. Copy and paste the XYZT and Frame columns from the TrackMate spots CSV file into the corresponding orange columns in the spreadsheet. Enter the numerical aperture (NA) and emission wavelength used during the experiment. Calculate the relative displacement in X, Y, and Z using the formula: Relative Displacement = Position - Position_Initial. Finally, plot the relative displacement over time to visualize the system's drift.

Identify visually the time when the displacement is lower than the resolution of the system. On this instrument it takes 120 min to reach its stability. Calculate the velocity, Velocity = (Displacement₂-Displacement₁)/T₂-T₁) and plot the velocity over time.

Calculate the average velocity before and after stabilisation and report the results in a table.

The warmup time for this instrument is approximately 2 hours. After the warmup period, the average displacement velocity is 14 nm/min, which falls within an acceptable range.

1. Place 4 µm diameter fluorescent beads (TetraSpec Fluorescent Microspheres Size Kit, mounted on a slide) on the microscope stage.

2. Center an isolated bead under a high-NA dry objective.

3. Crop the acquisition area to only visualize one bead but keep it large enough to anticipate a potential drift along the X and Y axes (100 um FOV should be enough)
   

4. Select an imaging channel appropriate for the fluorescent beads

   I typically use the Cy5 channel which is very bright and resistant to bleaching, yet this channel has a lower resolution.

5. Acquire a Z-stack at positions separated by 0 um (Drift), 1000 um and 10 000 um in both X and Y direction.

6. Repeat the acquisition 20 times

7. Acquire 3 datasets for each condition

1. Open your image in FIJI
2. If necessary, crop the image to focus on a single bead for better visualization.
3. Use the TrackMate plugin included in FIJI to detect and track spots over time Plugins\Tracking\TrackMate.
4. Apply Difference of Gaussians (DoG) spot detection with a detection size of 4 µm
5. Enable sub-pixel localization for increased accuracy
6. Click Preview to visualize the spot detection
7. Set a quality threshold (click and slide) high enough to detect a single spot per frame
8. Click Next and follow the detection and tracking process
9. Save the detected Spots coordinates as a CSV file for further analysis

Open the spreadsheet template XYZ_Repositioning Dispersion_Template.xlsx and fill in the orange cells. Copy and paste the XYZT and Frame columns from the TrackMate spots CSV file into the corresponding orange columns in the spreadsheet. Enter the numerical aperture (NA) and emission wavelength used during the experiment. Calculate the relative position in X, Y, and Z using the formula: Relative Position_Relative = Position - Position_Initial. Finally, plot the relative position over time to visualize the system's stage repositioning dispersion.

We observe an initial movement in X and Y that stabilises. Calculate the displacement 2D Displacement = Sqrt( (X₂-X₁)² + (Y₂-Y₁)²) ) and plot the 2D displacement over time. Calculate the resolution of your imaging configuration, Lateral Resolution = Lambda_Emission / 2*NA and plot the resolution over time (constant).

TThis experiment shows a significant initial displacement between Frame 0 and Frame 1, ranging from 1000 nm to 400 nm, which decreases to 70 nm by Frame 2. To quantify this variation, calculate the dispersion for each parameter using the formula: Dispersion = Standard Deviation (Value). Report the resultsin a table.

The system has a stage dispersion of 230 nm. The results are higher than expected because of the initial shift in the first frame that eventually stabilizes. Excluding the first frame significantly improves the measurements, reducing the repositioning dispersion to 40 nm. Further investigation is required to understand the underlying cause.

We observed a significant shift in the first frame, which was unexpected and invites further investigation. These variables can affect repositioning dispersion: i) Traveled distance, ii) Speed, iii) Acceleration, iv) Time, and v) Environment. We decided to test the first three.

Methodology

import os
import sys
import csv
from ij import IJ, WindowManager, Prefs
from ij.gui import GenericDialog
from java.io import File
from javax.swing import JOptionPane, JFileChooser  # Import JOptionPane and JFileChooser for dialogs
from fiji.plugin.trackmate import Model, Settings, TrackMate, SelectionModel, Logger
from fiji.plugin.trackmate.detection import DogDetectorFactory
from fiji.plugin.trackmate.tracking.jaqaman import SparseLAPTrackerFactory
from fiji.plugin.trackmate.features import FeatureFilter
from fiji.plugin.trackmate.features.track import TrackIndexAnalyzer
from fiji.plugin.trackmate.gui.displaysettings import DisplaySettingsIO
from fiji.plugin.trackmate.gui.displaysettings.DisplaySettings import TrackMateObject
from fiji.plugin.trackmate.visualization.table import TrackTableView
from fiji.plugin.trackmate.visualization.hyperstack import HyperStackDisplayer
from loci.plugins import BF
from loci.plugins.in import ImporterOptions

# Ensure UTF-8 encoding
reload(sys)
sys.setdefaultencoding('utf-8')

# Default settings in case there's no previous config file
TrackMateSettingsDefault = {
    'subpixel_localization': True,
    'spot_diameter': 4.0,  # Default 0.5 microns
    'threshold_value': 20.904,  # Default threshold value
    'apply_median_filtering': False
}


# Define keys for preferences
PREFS_KEY_SUBPIXEL = "BatchTrackmate.subpixel_localization"
PREFS_KEY_DIAMETER = "BatchTrackmate.spot_diameter"
PREFS_KEY_THRESHOLD = "BatchTrackmate.threshold_value"
PREFS_KEY_MEDIAN = "BatchTrackmate.apply_median_filtering"
PREFS_KEY_PROCESS_ALL_IMAGES = "BatchTrackmate.process_all_images"

def read_settings_from_prefs():
    subpixel_localization_str = Prefs.get(PREFS_KEY_SUBPIXEL, "True")
    spot_diameter_str = Prefs.get(PREFS_KEY_DIAMETER, "4.0")
    threshold_value_str = Prefs.get(PREFS_KEY_THRESHOLD, "20.904")
    apply_median_filtering_str = Prefs.get(PREFS_KEY_MEDIAN, "False")
    process_all_images_str = Prefs.get(PREFS_KEY_PROCESS_ALL_IMAGES, "False")

    TrackMateSettingsStored = {
        'subpixel_localization': subpixel_localization_str.lower() == "true",
        'spot_diameter': float(spot_diameter_str),
        'threshold_value': float(threshold_value_str),
        'apply_median_filtering': apply_median_filtering_str.lower() == "true",
        'process_all_images': process_all_images_str.lower() == "true"
    }
    return TrackMateSettingsStored

def save_settings_to_prefs(TrackMateSettings):
    Prefs.set(PREFS_KEY_SUBPIXEL, str(TrackMateSettings['subpixel_localization']))
    Prefs.set(PREFS_KEY_DIAMETER, str(TrackMateSettings['spot_diameter']))
    Prefs.set(PREFS_KEY_THRESHOLD, str(TrackMateSettings['threshold_value']))
    Prefs.set(PREFS_KEY_MEDIAN, str(TrackMateSettings['apply_median_filtering']))
    Prefs.set(PREFS_KEY_PROCESS_ALL_IMAGES, str(TrackMateSettings['process_all_images']))
    Prefs.savePreferences()

def display_dialog(TrackMateSettings, imp, NbDetectedSpots):
    NbTimepoints = imp.getNFrames()
    Max_Intensity = imp.getProcessor().getMax()
    gd = GenericDialog("Batch TrackMate Parameters")

    # Add fields to the dialog, prefill with last saved values
    gd.addMessage("=== Detection Settings ===")
    gd.addCheckbox("Enable Subpixel Localization", TrackMateSettings['subpixel_localization'])
    gd.addNumericField("Spot Diameter (microns):", TrackMateSettings['spot_diameter'], 2, 15, "")
    gd.addSlider("Threshold Value:", 0, Max_Intensity, TrackMateSettings['threshold_value'])
    gd.addCheckbox("Apply Median Filtering", TrackMateSettings['apply_median_filtering'])
    gd.addCheckbox("Process All Images", TrackMateSettings['process_all_images'])
    if NbDetectedSpots == False:
        NbDetectedSpots = 0
    if NbDetectedSpots > NbTimepoints:
        gd.addMessage("Nb of detected spots too high. Increase the threshold.")
    elif NbDetectedSpots < NbTimepoints:
        gd.addMessage("Nb of detected spots too low. Decrease the threshold.")
    elif NbDetectedSpots == NbTimepoints:
        gd.addMessage("The nb of detected spots matches the nb of timepoints.")
    gd.addMessage("Nb of detected spots: " + str(NbDetectedSpots) + ". Nb of Timepoints: " + str(NbTimepoints) + ".")
    gd.addButton("Test Detection", lambda event: Test_Detection(event, gd))
    # Show the dialog
    gd.showDialog()
    User_Click = None  # Initialize User_Click within the function
    if gd.wasOKed():
        User_Click = "OK"
    elif gd.wasCanceled():
        User_Click = "Cancel"

    # Get user inputs from the dialog
    TrackMateSettingsUser = {}
    TrackMateSettingsUser['subpixel_localization'] = gd.getNextBoolean()
    TrackMateSettingsUser['spot_diameter'] = gd.getNextNumber()
    TrackMateSettingsUser['threshold_value'] = gd.getNextNumber()
    TrackMateSettingsUser['apply_median_filtering'] = gd.getNextBoolean()
    TrackMateSettingsUser['process_all_images'] = gd.getNextBoolean()
    
    return TrackMateSettingsUser, User_Click, NbDetectedSpots, NbTimepoints

def run_trackmate_detection(imp, TrackMateSettings, SaveCSV):
    model = Model()
    model.setLogger(Logger.IJ_LOGGER)  # Send messages to the ImageJ log window

    # Configure TrackMate for this image
    trackmate_settings = Settings(imp)
    trackmate_settings.detectorFactory = DogDetectorFactory()
    trackmate_settings.detectorSettings = {
        'DO_SUBPIXEL_LOCALIZATION': TrackMateSettings['subpixel_localization'],
        'RADIUS': TrackMateSettings['spot_diameter'] / 2,
        'TARGET_CHANNEL': 1,
        'THRESHOLD': TrackMateSettings['threshold_value'],
        'DO_MEDIAN_FILTERING': TrackMateSettings['apply_median_filtering'],
    }

    # Configure tracker
    trackmate_settings.trackerFactory = SparseLAPTrackerFactory()
    trackmate_settings.trackerSettings = trackmate_settings.trackerFactory.getDefaultSettings()
    trackmate_settings.trackerSettings['ALLOW_TRACK_SPLITTING'] = True
    trackmate_settings.trackerSettings['ALLOW_TRACK_MERGING'] = True

    # Add all known feature analyzers to compute track statistics
    trackmate_settings.addAllAnalyzers()

    # Configure track filters
    track_filter = FeatureFilter('TRACK_DISPLACEMENT', 10, False)
    trackmate_settings.addTrackFilter(track_filter)

    # Instantiate plugin
    trackmate = TrackMate(model, trackmate_settings)

    # Process
    if not trackmate.checkInput():
        IJ.log("Error checking input: " + trackmate.getErrorMessage())
        return False
    if not trackmate.process():
        IJ.log("Error during processing: " + trackmate.getErrorMessage())
        return False

    # Display results
    selection_model = SelectionModel(model)

    ds = DisplaySettingsIO.readUserDefault()
    ds.setTrackColorBy(TrackMateObject.TRACKS, TrackIndexAnalyzer.TRACK_INDEX)
    ds.setSpotColorBy(TrackMateObject.TRACKS, TrackIndexAnalyzer.TRACK_INDEX)

    displayer = HyperStackDisplayer(model, selection_model, imp, ds)
    displayer.render()
    displayer.refresh()

    # Display the number of detected spots and the number of timepoints
    NbDetectedSpots = model.getSpots().getNSpots(True)
    if SaveCSV:
        # Export spot table as CSV
        spot_table = TrackTableView.createSpotTable(model, ds)
        filename = imp.getTitle()
        output_csv_path = generate_unique_filepath(output_dir, filename, ".csv")
        spot_table.exportToCsv(File(output_csv_path))
    
    return NbDetectedSpots

def generate_unique_filepath(directory, base_name, extension):
    """
    Generates a unique filename by appending a numeric suffix.
    
    :param directory: The directory where the file will be saved.
    :param base_name: The base name of the file (without the suffix and extension).
    :param extension: The file extension (e.g., '.txt').
    :return: A unique filename with an incremented suffix.
    """
    suffix = 1
    while True:
        filename = "{}_Spots-{:03d}{}".format(base_name, suffix, extension)
        full_path = os.path.join(directory, filename)
        if not os.path.exists(full_path):
            return full_path
        suffix += 1

User_Click = None  # Initialize User_Click
processed_files = []  # List to keep track of processed files

def Test_Detection(event, gd):
    global User_Click
    User_Click = "Retry"
    gd.dispose()  # Close the dialog

def select_folder():
    """
    Prompts the user to select a folder using JFileChooser.
    
    :return: The selected folder path or None if no folder was selected.
    """
    desktop_dir = os.path.join(os.path.expanduser("~"), "Desktop")
    chooser = JFileChooser(desktop_dir)
    chooser.setFileSelectionMode(JFileChooser.DIRECTORIES_ONLY)
    chooser.setDialogTitle("Choose a directory containing the images to process")
    return_value = chooser.showOpenDialog(None)
    if return_value == JFileChooser.APPROVE_OPTION:
        return chooser.getSelectedFile().getAbsolutePath()
    return None

def open_image_with_bioformats(file_path):
    """
    Opens an image using Bio-Formats and returns the ImagePlus object.
    
    :param file_path: Path to the image file.
    :return: ImagePlus object or None if the file could not be opened.
    """
    options = ImporterOptions()
    options.setId(file_path)
    try:
        imps = BF.openImagePlus(options)
        if imps is not None and len(imps) > 0:
            return imps[0]
        else:
            IJ.log("No images found in the file: " + file_path)
            return None
    except Exception as e:
        IJ.log("Error opening file with Bio-Formats: " + str(e))
        return None

# Main processing loop
folder_path = None
image_from_folder = False

while True:
    TrackMateSettingsStored = read_settings_from_prefs()
    imp = WindowManager.getCurrentImage()
    
    # If no image is open, prompt the user to select a folder
    if imp is None:
        folder_path = select_folder()
        if folder_path is not None:
            image_files = [f for f in os.listdir(folder_path) if f.lower().endswith(('.tif', '.tiff', '.jpg', '.jpeg', '.png', '.czi', '.nd2', '.lif', '.lsm', '.ome.tif', '.ome.tiff'))]
            if not image_files:
                JOptionPane.showMessageDialog(None, "No valid image files found in the selected folder.", "Error", JOptionPane.ERROR_MESSAGE)
                sys.exit(1)
            imp = open_image_with_bioformats(os.path.join(folder_path, image_files[0]))
            if imp is None:
                JOptionPane.showMessageDialog(None, "Failed to open the first image in the folder.", "Error", JOptionPane.ERROR_MESSAGE)
                sys.exit(1)
            image_from_folder = True
        else:
            JOptionPane.showMessageDialog(None, "No folder selected. Exiting.", "Error", JOptionPane.ERROR_MESSAGE)
            sys.exit(1)

    NbDetectedSpots = run_trackmate_detection(imp, TrackMateSettingsStored, SaveCSV=False)
    TrackMateSettingsUser, User_Click, NbDetectedSpots, NbTimepoints = display_dialog(TrackMateSettingsStored, imp, NbDetectedSpots)
    save_settings_to_prefs(TrackMateSettingsUser)
    
    if User_Click == "OK" and NbDetectedSpots == NbTimepoints and TrackMateSettingsStored == TrackMateSettingsUser:
        break
    elif User_Click == "Cancel":
        IJ.log("User canceled the operation.")
        if image_from_folder:
            imp.close()
        sys.exit(1)


# Create output directory on the user's desktop
output_dir = os.path.join(os.path.expanduser("~"), "Desktop", "Output")
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# Initialize a list to store data for each processed file
processed_data = []


if TrackMateSettingsUser['process_all_images']:
    if folder_path is not None:
        for image_file in image_files:
            imp.close()
            imp = open_image_with_bioformats(os.path.join(folder_path, image_file))
            if imp is not None:
               NbDetectedSpots = run_trackmate_detection(imp, TrackMateSettingsStored, SaveCSV=True)
               processed_files.append(image_file)
               processed_data.append({
                    "Filename": image_file,
                    "NbDetectedSpots": NbDetectedSpots,
                    "NbTimepoints": imp.getNFrames()
                    })
               imp.close()
    else:
        image_ids = WindowManager.getIDList()
        if image_ids is not None:
            for image_id in image_ids:
                imp = WindowManager.getImage(image_id)
                NbDetectedSpots = run_trackmate_detection(imp, TrackMateSettingsStored, SaveCSV=True)
                processed_files.append(imp.getTitle())
                processed_data.append({
                    "Filename": imp.getTitle(),
                    "NbDetectedSpots": NbDetectedSpots,
                    "NbTimepoints": imp.getNFrames()
                })
else:
    NbDetectedSpots = run_trackmate_detection(imp, TrackMateSettingsStored, SaveCSV=True)
    processed_files.append(imp.getTitle())
    processed_data.append({
        "Filename": imp.getTitle(),
        "NbDetectedSpots": NbDetectedSpots,
        "NbTimepoints": imp.getNFrames()
    })

# Display a message dialog at the end of processing
message = "Processing Complete!\nNumber of files processed: {}\nOutput files are saved in: {}".format(len(processed_files), output_dir)
JOptionPane.showMessageDialog(None, message, "Processing Complete", JOptionPane.INFORMATION_MESSAGE)

# Close the first image if it was opened from a folder
if image_from_folder:
    imp.close()
    
# Display summary in the log at the end
IJ.log("=== Summary of Processed Files ===")
for data in processed_data:
    IJ.log("Filename: {} | Detected spots: {} | Number of timepoints: {}".format(
        data["Filename"], data["NbDetectedSpots"], data["NbTimepoints"]
    ))

# Load and install necessary libraries at the beginning
rm(list = ls())

# Check and install required packages if they are not already installed
if (!require(dplyr)) install.packages("dplyr", dependencies = TRUE)
if (!require(stringr)) install.packages("stringr", dependencies = TRUE)
if (!require(ggplot2)) install.packages("ggplot2", dependencies = TRUE)
if (!require(tcltk)) install.packages("tcltk", dependencies = TRUE)

# Load libraries
library(dplyr)
library(stringr)
library(ggplot2)
library(tcltk)

# Set default input and output directories
default_input_dir <- file.path(Sys.getenv("USERPROFILE"), "Desktop")  # Default to "Output" on Desktop
InputFolder <- tclvalue(tkchooseDirectory(initialdir = default_input_dir, title = "Select a folder containing CSV files"))

#InputFolder <- default_input_dir  # Use default folder

# Specify the Output folder (this is a fixed folder on the Desktop)
OutputFolder <- file.path(Sys.getenv("USERPROFILE"), "Desktop", "Output")
if (!dir.exists(OutputFolder)) dir.create(OutputFolder, recursive = TRUE)

# List all CSV files in the folder
csv_files <- list.files(path = InputFolder, pattern = "\\.csv$", full.names = TRUE)
if (length(csv_files) == 0) stop("No CSV files found in the selected directory.")

header <- names(read.csv(csv_files[1], nrows = 1))

# Function to clean the filename
clean_filename <- function(filename) {
  filename_parts <- strsplit(filename, " - ")[[1]]
  cleaned_filename <- sub("\\.czi", "", filename_parts[1])
  cleaned_filename <- sub("_spots", "", cleaned_filename)
  cleaned_filename <- sub("\\.csv", "", cleaned_filename)
  return(cleaned_filename)
}

# Read and merge all CSV files
merged_data <- csv_files %>%
  lapply(function(file) {
    data <- read.csv(file, skip = 4, header = FALSE)
    colnames(data) <- header
    data <- data %>% arrange(FRAME)
    
    # Clean and add source file info
    filename <- basename(file)
    data$SourceFile <- clean_filename(filename)
    
    # Extract variables from the filename
    filename_parts <- strsplit(clean_filename(filename), "_")[[1]]
    for (i in seq_along(filename_parts)) {
      variable_name <- paste("Variable-", sprintf("%03d", i), sep = "")
      data[[variable_name]] <- filename_parts[i]
    }
    
    # Add time columns if available
    if ("POSITION_T" %in% colnames(data)) {
      data$`Time (sec)` <- round(data$POSITION_T, 0)
      data$`Time (min)` <- round(data$`Time (sec)` / 60, 2)
    }
    
    # Calculate displacement columns (X, Y, Z, 3D, 2D)
    if ("POSITION_X" %in% colnames(data)) {
      first_value <- data$POSITION_X[data$FRAME == 0][[1]]
      data$`X (nm)` <- (data$POSITION_X - first_value) * 1000
    }
    if ("POSITION_Y" %in% colnames(data)) {
      first_value <- data$POSITION_Y[data$FRAME == 0][[1]]
      data$`Y (nm)` <- (data$POSITION_Y - first_value) * 1000
    }
    if ("POSITION_Z" %in% colnames(data)) {
      first_value <- data$POSITION_Z[data$FRAME == 0][[1]]
      data$`Z (nm)` <- (data$POSITION_Z - first_value) * 1000
    }
    
    # Calculate displacement (3D and 2D)
    if (all(c("X (nm)", "Y (nm)", "Z (nm)") %in% colnames(data))) {
      data$`Displacement 3D (nm)` <- sqrt(
        diff(c(0, data$`X (nm)`))^2 + 
          diff(c(0, data$`Y (nm)`))^2 + 
          diff(c(0, data$`Z (nm)`))^2
      )
      data$`Displacement 3D (nm)`[1] <- 0
    }
    if (all(c("X (nm)", "Y (nm)") %in% colnames(data))) {
      data$`Displacement 2D (nm)` <- sqrt(
        diff(c(0, data$`X (nm)`))^2 + 
          diff(c(0, data$`Y (nm)`))^2
      )
      data$`Displacement 2D (nm)`[1] <- 0
    }
    
    return(data)
  }) %>%
  bind_rows()

# Save the merged data to a new CSV file
output_file <- "XYZ_Repositioning_Merged-Data.csv"
write.csv(merged_data, file = file.path(OutputFolder, output_file), row.names = FALSE)

cat("All CSV files have been merged and saved to", output_file, "\n")

 Repositioning Dispersion: Impact of Traveled Distance

Results

Conclusion

Repositioning Dispersion: Impact of Speed and Acceleration

Results

Conclusion

What about the initial shift ?

Right, I almost forgot about that. See below.

Results

Conclusion

1. Place 4 µm diameter fluorescent beads (TetraSpec Fluorescent Microspheres Size Kit, mounted on a slide) on the microscope stage.

2. Center an isolated bead under a high-NA dry objective.

3. Crop the acquisition area to only visualize one bead but keep it large enough to anticipate a potential chromatic shift along the XY and Z axes
   

4. Acquire a multi-channel Z-stack
   

   I usually acquire all available channels, with the understanding that no more than seven channels can be processed simultaneously. If you have more than 7 channels you can either split them into sets smaller or equal to 7 keeping one reference channel in all sets.
   

5. Acquire 3 datasets
   

6. Repeat for each objective

You should have acquired several multi-channel images that now need processing to yield meaningful results. To process them, use the Channel Co-registration analysis feature of the MetroloJ_QC plugin for FIJI. For more information about the MetroloJ_QC plugin please refer to manual available at the MontpellierRessourcesImagerie on GitHub.

1. Open FIJI.
2. Load your image by dragging it into the FIJI bar.
3. Launch the MetroloJ_QC plugin by navigating to Plugins > MetroloJ QC.
4. Click on Channel Co-registration Report.
5. Enter a Title for your report.
6. Type in your Name.
7. Click Microscope Acquisition Parameters and enter:
   1. The imaging modality (Widefield, Confocal, Spinning-Disk, Multi-photon) 
   2. The Numerical Aperture of the Objective used.
   3. The Refractive index of the imersion media. (Air: 1.0; Water 1.33; Oil: 1.515)
      
   4. The Emission Wavelength in nm for each channel.

8. Click OK

9. Click Bead Detection Options
   1. Bead Detection Threshold: Legacy
   2. Center Method: Legacy Fit ellipses

10. Enable Apply Tolerance to the Report and Reject co-registration ratio above 1.

11. Click File Save Options.

12. Select Save Results as PDF Reports.

13. Select Save Results as spreadsheets.

14. Select Save result images.
15. Select Open individual pdf report(s).

16. Click OK.

17. Repeat steps 4 through 16 for each image you have acquired

This will generate detailed results stored in a folder named Processed. The Processed folder will be located in the same directory as the original images, with each report and result saved in its own sub-folder.

Here is a script to process channel co-alignment results generated by MetroloJ_QC Process Ch Co-registration Results.R.

To use it, simply drag and drop the file into the R interface. You may also open it with RStudio and Click the Source button. The script will:

• Prompt the user to select an input folder.
• Load all _results.xls files located in the selected folder and subfolder
  
• Process and merge all results generated by the MetroloJ_QC Channel Co-registration.
• Split the filenames using the _ character as a separator and create columns named Variable-001, Variable-002, etc. This will help in organizing the data if the filenames are formatted as Variable-A_Variable-B, for example.
• Save the result as Channel_Co-registration_Merged-Data.csv on an Output folder on the user's Desktop

# Clear the workspace
rm(list = ls())

# Ensure required packages are installed and loaded
if (!require(tcltk)) install.packages("tcltk", dependencies = TRUE)
library(tcltk)
if (!require(reshape2)) install.packages("reshape2", dependencies = TRUE)
library(reshape2)


# Clear the workspace
rm(list = ls())

# Set default input and output directories
default_input_dir <- file.path(Sys.getenv("USERPROFILE"), "Desktop", "Input")  # Default to "Input" on Desktop
InputFolder <- default_input_dir  # Use default folder
InputFolder <- tclvalue(tkchooseDirectory(initialdir = default_input_dir, title = "Select a folder containing XLS files"))

OutputFolder <- file.path(Sys.getenv("USERPROFILE"), "Desktop", "Output")
if (!dir.exists(OutputFolder)) dir.create(OutputFolder, recursive = TRUE)

# List all XLS files with "_results.xls" in the folder and subfolders
xls_files <- list.files(path = InputFolder, pattern = "_results.xls$", full.names = TRUE, recursive = TRUE)
if (length(xls_files) == 0) stop("No _results.xls files found in the selected directory.")

# Initialize lists for storing data
ratio_data_tables <- list()
pixel_shift_tables <- list()

# Process each XLS file
for (file in xls_files) {
  lines <- readLines(file, encoding = "latin1")
  
  # Extract "Ratios" data
  ratios_line_index <- grep("Ratios", lines)
  if (length(ratios_line_index) == 0) {
    cat("No 'Ratios' line found in file:", file, "\n")
    next
  }
  
  start_row <- ratios_line_index[1] + 1
  relevant_lines <- lines[start_row:length(lines)]
  
  channel_start_index <- which(grepl("^Channel", relevant_lines))[1]
  if (is.na(channel_start_index)) {
    cat("No rows starting with 'Channel' found after 'Ratios' in file:", file, "\n")
    next
  }
  
  contiguous_rows <- relevant_lines[channel_start_index]
  for (i in (channel_start_index + 1):length(relevant_lines)) {
    if (grepl("^Channel", relevant_lines[i])) {
      contiguous_rows <- c(contiguous_rows, relevant_lines[i])
    } else {
      break
    }
  }
  
  if (length(contiguous_rows) > 0) {
    # Read the ratio_data table from contiguous rows
    ratio_data <- read.table(text = paste(contiguous_rows, collapse = "\n"), sep = "\t", header = FALSE, stringsAsFactors = FALSE)
    
    # Extract channel names and matrix of values
    channels <- ratio_data$V1
    values <- ratio_data[, -1]
    colnames(values) <- channels
    rownames(values) <- channels
    
    # Melt the matrix into a long-format table
    ratio_data_formatted <- melt(as.matrix(values), varnames = c("Channel_1", "Channel_2"), value.name = "Ratio")
    ratio_data_formatted$Pairing <- paste(ratio_data_formatted$Channel_1, "x", ratio_data_formatted$Channel_2)
    
    # Add filename and split into parts
    filename_no_extension <- tools::file_path_sans_ext(basename(file))
    ratio_data_formatted$Filename <- filename_no_extension
    
    filename_parts <- strsplit(ratio_data_formatted$Filename, "_")
    max_parts <- max(sapply(filename_parts, length))
    for (i in 1:max_parts) {
      new_column_name <- paste0("Variable-", sprintf("%03d", i))
      ratio_data_formatted[[new_column_name]] <- sapply(filename_parts, function(x) ifelse(length(x) >= i, x[i], NA))
    }
    
    # Add channel name pairings
    channel_names <- unlist(strsplit(ratio_data_formatted$`Variable-005`[1], "-"))
    ratio_data_formatted$Channel_Name_Pairing <- apply(ratio_data_formatted, 1, function(row) {
      ch1_index <- as.integer(gsub("Channel ", "", row["Channel_1"])) + 1
      ch2_index <- as.integer(gsub("Channel ", "", row["Channel_2"])) + 1
      paste(channel_names[ch1_index], "x", channel_names[ch2_index])
    })
    # Filter columns for ratio_data: Channel_1, Channel_2, Channel_Pairing, Ratio, Filename, Variable-001, Variable-002 ..., Channel_Name_Pairing
    ratio_data_formatted <- ratio_data_formatted[, c("Channel_1", "Channel_2", "Pairing", "Ratio", "Filename", 
                                                     paste0("Variable-", sprintf("%03d", 1:7)), "Channel_Name_Pairing")]
    
    ratio_data_tables[[basename(file)]] <- ratio_data_formatted
  } else {
    cat("No contiguous 'Channel' rows found to convert into a data frame.\n")
  }
  
  # Extract "Pixel Shifts" data
  pixelShifts_line_index <- grep("pixelShifts", lines)
  if (length(pixelShifts_line_index) > 0) {
    start_row_pixel_shifts <- pixelShifts_line_index[1] + 1
    relevant_pixel_shifts_lines <- lines[start_row_pixel_shifts:length(lines)]
    pixel_shift_rows <- character(0)
    
    for (line in relevant_pixel_shifts_lines) {
      columns <- strsplit(line, "\t")[[1]]
      if (length(columns) >= 2 && grepl("shift$", columns[2])) {
        pixel_shift_rows <- c(pixel_shift_rows, line)
      }
    }
    
    if (length(pixel_shift_rows) > 0) {
      pixel_shift_data <- read.table(text = paste(pixel_shift_rows, collapse = "\n"), sep = "\t", header = FALSE, stringsAsFactors = FALSE)
      # Set the first column name to "Channel"
      colnames(pixel_shift_data)[1] <- "Channel"
      # Extract unique channel names from the `Channel` column
      # Filter out empty strings and get unique channel values
      unique_channels <- unique(pixel_shift_data$Channel[pixel_shift_data$Channel != ""])
      
      # Update column names: Keep "Channel" and "Axis", then use `unique_channels` for the rest
      new_colnames <- c("Channel", "Axis", unique_channels)
      
      # Update the column names of the data frame
      colnames(pixel_shift_data) <- new_colnames
      
      for (i in 2:nrow(pixel_shift_data)) {
        if (is.na(pixel_shift_data$Channel[i]) || pixel_shift_data$Channel[i] == "") {
          pixel_shift_data$Channel[i] <- pixel_shift_data$Channel[i - 1]
        }
      }
      
      # Melt the pixel_shift_data into long format
      pixel_shift_long <- melt(pixel_shift_data, id.vars = c("Channel", "Axis"), variable.name = "Channel_2", value.name = "Pixel_Shift")
      colnames(pixel_shift_long)[1]<-"Channel_1"
      
      # Create the Channel Pairing column
      pixel_shift_long$Pairing <- paste(pixel_shift_long$Channel_1, "x", pixel_shift_long$Channel_2)
      
      # Filter data into separate X, Y, Z based on Axis
      x_shifts <- subset(pixel_shift_long, Axis == "X shift", select = c(Pairing, Pixel_Shift))
      y_shifts <- subset(pixel_shift_long, Axis == "Y shift", select = c(Pairing, Pixel_Shift))
      z_shifts <- subset(pixel_shift_long, Axis == "Z shift", select = c(Pairing, Pixel_Shift))
      
      # Rename Pixel_Shift columns for clarity
      colnames(x_shifts)[2] <- "X"
      colnames(y_shifts)[2] <- "Y"
      colnames(z_shifts)[2] <- "Z"
      
      # Merge X, Y, Z shifts into a single data frame
      pixel_shift_data_formatted <- Reduce(function(x, y) merge(x, y, by = "Pairing", all = TRUE),
                                           list(x_shifts, y_shifts, z_shifts))
      
      
      pixel_shift_data_formatted$Filename <- filename_no_extension
      for (i in 1:max_parts) {
        new_column_name <- paste0("Variable-", sprintf("%03d", i))
        pixel_shift_data_formatted[[new_column_name]] <- sapply(filename_parts, function(x) ifelse(length(x) >= i, x[i], NA))
      }
      
      
      # Split Variable-005 into individual channel names
      channel_names <- strsplit(pixel_shift_data_formatted$`Variable-005`, split = "-")
      
      # Create a function to match channel names to Channel_Pairing
      get_channel_pairing_name <- function(pairing, names_list) {
        # Extract the channel numbers from the pairing
        channel_indices <- as.numeric(gsub("Channel ", "", unlist(strsplit(pairing, " x "))))
        # Return the corresponding names from names_list
        paste(names_list[channel_indices + 1], collapse = " x ")
      }
      
      # Apply the function to generate the Channel_Pairing_Name column
      pixel_shift_data_formatted$Channel_Name_Pairing <- mapply(
        get_channel_pairing_name,
        pairing = pixel_shift_data_formatted$Pairing,
        MoreArgs = list(names_list = unlist(channel_names[1]))
      )
      # Split Pairing into Channel_1 and Channel_2 properly
      Single_channel_nb <- strsplit(pixel_shift_data_formatted$Pairing, split = " x ")
      
      # Assign Channel_1 and Channel_2
      pixel_shift_data_formatted$Channel_1 <- sapply(Single_channel_nb, `[`, 1)
      pixel_shift_data_formatted$Channel_2 <- sapply(Single_channel_nb, `[`, 2)
      
      
      
      pixel_shift_tables[[basename(file)]] <- pixel_shift_data_formatted
    }
  }
}

# Combine and save data
combined_ratio_data <- do.call(rbind, ratio_data_tables)
output_ratio_file <- file.path(OutputFolder, "Ratio_Data.csv")
#write.csv(combined_ratio_data, output_ratio_file, row.names = FALSE)

combined_pixel_shift_data <- do.call(rbind, pixel_shift_tables)
output_pixel_shift_file <- file.path(OutputFolder, "Pixel_Shift_Data.csv")
#write.csv(combined_pixel_shift_data, output_pixel_shift_file, row.names = FALSE)

merge_columns <- c("Filename", "Channel_Name_Pairing", "Channel_1","Channel_2","Pairing", paste0("Variable-", sprintf("%03d", 1:7)))
combined_data <- merge(combined_ratio_data, combined_pixel_shift_data, by = merge_columns, all = TRUE, suffixes = c("_Ratio", "_PixelShift"))

# Split Pairing into Channel_1 and Channel_2 properly
Single_channel_Name <- strsplit(combined_data$Channel_Name_Pairing, split = " x ")

# Assign Channel_1 and Channel_2
combined_data$Channel_1_Name <- sapply(Single_channel_Name, `[`, 1)
combined_data$Channel_2_Name <- sapply(Single_channel_Name, `[`, 2)

# Function to create unique filenames by incrementing the suffix
get_unique_filename <- function(base_path) {
  count <- 0
  new_file <- base_path
  while (file.exists(new_file)) {
    count <- count + 1
    new_file <- paste0(sub("\\.csv$", "", base_path), sprintf("-%03d.csv", count))
  }
  return(new_file)
}

output_combined_file <- file.path(OutputFolder, "Channel_Co-registration_Merged-Data.csv")
output_combined_file <- get_unique_filename(output_combined_file)  # Ensure unique filename
write.csv(combined_data, output_combined_file, row.names = FALSE)

The following spreadsheet provides a dataset that can be manipulated with a pivot table to generate informative graphs and statistics Channel_Co-registration_Data_Template.xlsx.

 Plot the uniformity and centering accuracy for each objective.

Channel_Registration_With TrackMate.py

 Plot the uniformity and centering accuracy for each filter set. 

This specific instrument has a quad-band filter as well as individual filter cubes. We can plot the uniformity and centering accuracy per filter types. 

Data has several variables: Objective, filters.

looking at Objective The ratio are below 1 all is good

Looking at indivudal filters we see that individaul filter cubes have a higher ratio than the quad band

Here is an example of a produced spreadsheet Channel_Co-registration_Template.xlsx

Report the Pixel Shift Table For each objective and each filter combination. This table can be used to correct a multi-channel image by displacing the Channel 2 relative to the Channel 1 by the XYZ pixel coordinates.

 Cube Pixel Shift table