-
Notifications
You must be signed in to change notification settings - Fork 1
8. Image and Signal Analysis
- suite2p: Calcium imaging analysis.
- EZcalcium: Calcium imaging analysis. Motion correction, ROI detection, ROI refinement, and ROI matching.
- PatchWarp: Calcium imaging analysis. Corrects complex image distortions that slowly occur during long imaging sessions.
- Cellpose: Cell and nucleus segmentation
- ilastik: Leverage machine learning algorithms to easily segment, classify, track and count your cells or other experimental data.
Some image processing software only work with tiff files. This scripts reads microscope images in their original formats(e.g. czi, nd2, lsm) and saves them as tiff.
// saveastiff_batch.ijm
// Author: Jane Ling (yan_to.ling@kcl.ac.uk)
// Last modified: 2025-02-18 15:23
//
// This macro saves images in fileDir as tif in saveDir.
// Note that image file names cannot contain " ". Use "_" or "-" instead.
//
// INPUTS
// - microscope images (e.g. czi, nd2, lsm)
//
// OUTPUTS
// If not splitting channels, ...
// - (original filename)_raw.tif
// If splitting channels, ...
// - C1-(original filename)_raw.tif
// - C2-(original filename)_raw.tif
// - C3-(original filename)_raw.tif
// MODIFY HERE: Path of the image files to be processed, searches all sub-directories.
fileDir = "/PATH/TO/PROJECT/originals/";
saveDir = "/PATH/TO/PROJECT/outputs/";
// MODIFY HERE: end of filename and input file type (file extension) e.g.".lsm"
endString = ".lsm";
// MODIFY HERE: whether to split channels and save as separate tif
splitChannel = true;
///////////////////////////////////////////////////
setBatchMode(true);
print("Saving tif...");
walkFilesSaveTif(fileDir, saveDir, splitChannel, endString);
print("Done:)");
////////////////////////////////////////////////////
// Local Functions //
////////////////////////////////////////////////////
// Find all files in subdirs, split channels and save tiff
function walkFilesSaveTif(inputDir, outputDir, splitChannel, endString) {
list = getFileList(inputDir);
// create output folder it it does not exist
if (!File.exists(outputDir)) {
File.makeDirectory(outputDir);
}
for (i=0; i<list.length; i++) {
// If folder
if (endsWith(list[i], File.separator))
walkFilesSaveTif( inputDir+list[i], outputDir+list[i], splitChannel, endString);
// If image file
else if (endsWith(list[i], endString)) {
print (list[i]);
// for saving all outputs in the same folder
// splitChannels(inputDir, outputDir, list[i]);
// for creating a separate folder for each image file
outputDir1 = outputDir + File.getNameWithoutExtension(list[i]) + File.separator;
saveTif(inputDir, outputDir1, list[i], splitChannel);
}
}
} // end of walkFilesSaveTif
// Split channels and save as tiff
function saveTif(inputDir, outputDir, file, splitChannel) {
open(inputDir + file);
// create output folder it it does not exist
if (!File.exists(outputDir)) {
File.makeDirectory(outputDir);
}
if (splitChannel) {
// total number of channels
Stack.getDimensions(width, height, channels, slices, frames);
Nc = channels;
run("Split Channels");
for (i=0; i<Nc; i++){
title = getTitle();
title = substring(title, 0, lastIndexOf(title, ".")); // get first part of file name
title = title + "_raw";
// if already saved, close
if (File.exists(outputDir + title + ".tif")) {
close();
}
else { // otherwise, save as tif
saveAs("Tiff", outputDir + title);
close();
}
}
}
else {
title = getTitle();
title = substring(title, 0, lastIndexOf(title, ".")); // get first part of file name
title = title + "_raw";
// if already saved, close
if (File.exists(outputDir + title + ".tif")) {
close();
}
else { // otherwise, save as tif
saveAs("Tiff", outputDir + title);
close();
}
}
} // end of saveTifThis script saves ROI boundaries as a list of coordinates in csv, which could be loaded later for processing.
// export_ROIs_as_csv.ijm
// Author: Jane Ling (yan_to.ling@kcl.ac.uk)
// Last modified: 2025-02-18 15:23
//
// This macro saves multiple ROIs from ROI Manager as csv
//
// Instructions
// 1. Open the image
// 2. Use any of the selection or line tool to draw an ROI.
// 3. Go to "Edit" -> "Selection" -> "Add to Manager".
// 4. Check that in "More >>" -> "Options", the checkbox "Associate "Show All" ROIs with slices" is checked.
// 5. Draw all other ROIs in the image/stack and add them to ROI Manager.
// 6. Go to"File" -> "Open".
// 7. Open the ImageJ macro file export_ROIs_as_csv.ijm
// 8. Modify fileDir to the saving directory on your computer.
// 9. Click Run to run the script. (Do NOT close the image before running the script.)
//
// INPUTS
// - ROIs from ROI Manager
//
// OUTPUTS
// - XY_(original filename)_1.csv
// - XY_(original filename)_2.csv
// - ...
// - (original filename)-ROIset.zip
print("Saving ROIs...");
// MODIFY HERE: Path to saving directory
fileDir = "/PATH/TO/PROJECT/outputs/";
// MODIFY HERE: Whether to create new folder under the saving directory with the name of the image
newFolder = true;
///////////////////////////////////////////////////
title = getTitle();
title = substring(title, 0, lastIndexOf(title, ".")); // get first part of file name
if (newFolder) {
fileDir = fileDir + title + File.separator;
// create output folder it it does not exist
if (!File.exists(fileDir)) {
File.makeDirectory(fileDir);
}
}
for (i=0; i<RoiManager.size; i++) {
roiManager("Select", i);
results_csv_path = fileDir + "XY_" + title + "_" + (i+1) + ".csv";
print(results_csv_path);
Roi.getCoordinates(x, y);
slice = getSliceNumber();
Stack.getPosition(channel, slice, frame);
run("Clear Results");
for (k=0; k<x.length; k++) {
setResult("X",k,x[k]);
setResult("Y",k,y[k]);
setResult("Z",k,slice);
}
updateResults();
selectWindow("Results");
saveAs("Results", results_csv_path);
roiManager("save", fileDir + title + "-ROIset.zip");
}
print("Done:)");Here are some materials for those who want to use MATLAB for image analysis.
- Self-paced online courses
- Analysis of microscopy images
- Medical Imaging Toolbox™
- Bioformats Image Toolbox
- Hausdorff distance
% close all figures
close all
% clear workspace
clear
% clear command window
clc
% Change the current folder to the folder of this m-file.
if(~isdeployed)
cd(fileparts(matlab.desktop.editor.getActiveFilename));
end
Bio-Formats is a convenient toolbox for reading biomedical images in various file formats. It has been installed on Clarke for Matlab under the folder /usr/local/bin/MATLAB/bfmatlab.
To load any data file, use the following lines:
% path to bioformat toolbox % change path if necessary
addpath('/PATH/TO/FOLDER/bfmatlab/')
ops.filename = "/PATH/TO/FOLDER/filename";
[im_data, ops] = loadBioFormats(ops);%% [im_data, ops] = loadBioFormats(ops, varargin)
% Loads data using bioformats toolbox and reshapes data.
%
% Author: Jane Ling (yan_to.ling@kcl.ac.uk)
% Last updated: 2025-02-14 12:40
%
% USAGE:
% 1) [im_data, ops] = loadBioFormats(ops);
% 2) [im_data, ops] = loadBioFormats(ops, verbose);
%
% INPUTS:
% - ops
% (struct) ops.filename defining the path to image file
% - verbose
% (bool) whether to pring information about dataset or not.
% (Default) false
%
% OUTPUTS:
% - im_data
% (numeric) 5D matrix if dimension order starts with 'XY', singular
% dimensions would be omitted.
% (cell array) interleaved frames otherwise
% - ops
% (struct) parameters of the stack
function [im_data, ops] = loadBioFormats(ops, varargin)
if nargin == 2
verbose = varargin{1};
else
verbose = 0;
end
% Bio-Formats Toolbox
[~,~] = bfCheckJavaPath(); % added such that path to Bio-Formats Toolbox is know
data = bfopen(ops.filename); % load data
im_data = data{1,1};
im_data = im_data(:,1); % image data
omeMeta = data{1,4}; % metadata
ops.imageIndex = 0;
% size of stack in each dimension
ops.Nx = omeMeta.getPixelsSizeX(ops.imageIndex).getValue(); % image width, pixels
ops.Ny = omeMeta.getPixelsSizeY(ops.imageIndex).getValue(); % image height, pixels
ops.Nz = omeMeta.getPixelsSizeZ(ops.imageIndex).getValue(); % number of Z slices
ops.Nc = omeMeta.getPixelsSizeC(ops.imageIndex).getValue(); % number of channels
ops.Nt = omeMeta.getPixelsSizeT(ops.imageIndex).getValue(); % number of time points
% dimension order in data
ops.dimOrder = char(omeMeta.getPixelsDimensionOrder(ops.imageIndex).getValue());
% size of pixels (voxels)
ops.sizeX = double(omeMeta.getPixelsPhysicalSizeX(ops.imageIndex).value(ome.units.UNITS.MICROMETER)); % [microns]
ops.sizeY = double(omeMeta.getPixelsPhysicalSizeY(ops.imageIndex).value(ome.units.UNITS.MICROMETER)); % [microns]
% axes
ops.x = (0:ops.Nx-1)*ops.sizeX;
ops.y = (0:ops.Ny-1)*ops.sizeY;
if verbose
fprintf('Size of image = %d x %d. \nNumber of Z slices = %d. \nNumber of channels = %d. \nNumber of time points = %d. \n', ops.Nx, ops.Ny, ops.Nz, ops.Nc, ops.Nt)
fprintf('Dimension order = %s \n', ops.dimOrder)
end
% for z-stack
if ops.Nz > 1
try
ops.sizeZ = double(omeMeta.getPixelsPhysicalSizeZ(ops.imageIndex).value(ome.units.UNITS.MICROMETER)); % [microns]
ops.z = (0:ops.Nz-1)*ops.sizeZ;
if verbose
fprintf('Voxel size = %.3f µm x %.3f µm x %.3f µm \n', ops.sizeX, ops.sizeY, ops.sizeZ)
end
catch
if verbose
fprintf('Pixel size = %.3f µm x %.3f µm \n', ops.sizeX, ops.sizeY)
end
warning('Could not read sizeZ from omeMeta.')
end
else
if verbose
fprintf('Pixel size = %.3f µm x %.3f µm \n', ops.sizeX, ops.sizeY)
end
end
% for time stack
if ops.Nt > 1
try
ops.sizeT = omeMeta.getPixelsTimeIncrement(ops.imageIndex).value(ome.units.UNITS.SECOND);
ops.t = (0:ops.Nt-1)*ops.sizeT;
if verbose
fprintf('Pixel time imcrement = %.3f s \n', ops.sizeT)
end
catch
warning('Dimension order may be incorrect.')
end
end
% reshaping
if strcmp(ops.dimOrder(1:2),'XY')
im_data = cell2mat(im_data);
im_data = reshape(im_data,ops.Ny,ops.Nc*ops.Nz*ops.Nt,ops.Nx);
im_data = permute(im_data, [1,3,2]);
N = [ops.Ny,ops.Nx, 0, 0, 0];
N(strfind(ops.dimOrder,'C')) = ops.Nc;
N(strfind(ops.dimOrder,'Z')) = ops.Nz;
N(strfind(ops.dimOrder,'T')) = ops.Nt;
im_data = reshape(im_data,N(1),N(2),N(3),N(4),N(5));
im_data = squeeze(im_data);
ops.dimOrder = ops.dimOrder(N~=1);
if verbose
fprintf('Dimension order of im_data after reshaping = %s \n', ops.dimOrder)
end
end
endIf you see and error that says Unrecognized function or variable, you may check if the folder has been added to path in HOME -> Set Path.
If you see an error related to java memory size, try to increase the java heap size following the steps describe in MATLAB Preferences (R2010a+)
For additional information, read Using Bio-Formats in MATLAB
| Function | Description |
|---|---|
| bfopen | Loads data using Bio-Formats toolbox |
| loadBioFormats | Custom function to load bioformats |
| load | Load variables from file (.mat) into workspace |
| imread | Read image from graphics file |
| readtable | Create table from file. For loading .xlsx or .csv |
| Function | Description |
|---|---|
| writetable | Write table to file |
| save | Save variables from workspace to file |
| saveas | Save figure to specific file format |
| imwrite | Write image to graphics file |
| saveastiff | Save multipage TIFF stack |
| Function | Description |
|---|---|
| fullfile | Build full file name from parts |
| strsplit | Split string or character vector at specified delimiter |
| strcat | Concatenate strings horizontally |
| filesep | File separator for current platform |
| fileparts | Get parts of file name |
| sprintf | Format data into string or character vector |
| Function | Description |
|---|---|
| double | Double-precision arrays |
| single | Single-precision arrays |
| int8 | 8-bit signed integer arrays |
| logical | Convert numeric values to logicals |
| string | String array |
| struct | Structure array |
| table | Table array with named variables that can contain different types |
| Function | Description |
|---|---|
| cell2mat | Convert cell array to ordinary array of the underlying data type |
| num2cell | Convert array to cell array with consistently sized cells |
| reshape | Reshape array by rearranging existing elements |
| permute | Permute array dimensions |
| repmat | Repeat copies of array |
| interp1 | 1-D data interpolation (table lookup). See also interp2, interp3, interpn |
| length | Length of largest array dimension |
| size | Array size |
| squeeze | Remove dimensions of length 1 |
| vertcat | Concatenate arrays vertically |
| horzcat | Concatenate arrays horizontally |
| Function | Description |
|---|---|
| medfilt2 | 2-D median filtering. See also medfilt3 |
| imgaussfilt | 2-D Gaussian filtering of images. See also imgaussfilt3 |
| fspecial | Create predefined 2-D filter |
| imfilter | N-D filtering of multidimensional images |
| imadjustn | Adjust intensity values in N-D volumetric image |
| Function | Description |
|---|---|
| filter | 1-D digital filter |
| movmean | Moving mean |
| movmedian | Moving median |
Zero padding helps to eliminate boundary/border effects of filtering by filling in the off-the-edge image pixels. Here is an example of zero padding by mirror-reflecting the array across the array border. For other common padding options, you may refer to the builtin imfilter function.
%% zero padding: mirror-reflected boundary
function data = zero_padding(data, zeropad_len)
Z1 = flipud(data(1:zeropad_len, :));
Z2 = flipud(data(end-zeropad_len+1:end,:));
data = [Z1; data; Z2];
end
function data = remv_zero_padding(data, zeropad_len)
data = data(zeropad_len+1:length(data)-zeropad_len,:);
end%% Sliding window filter
function baseline_pc_median = sliding_window_filter(data, baseline_percentage, window)
tic;
disp('Applying sliding window filter...')
zeropad_len = window;
data = zero_padding(data, zeropad_len);
T = size(data,1);
baseline_pc_median = zeros(size(data));
for k=1:T %for all timepoints
kymo_sample = data(max(1,k-window/2):min(T,k+window/2),:); %take window around each timepoint
sortedVals = sort(kymo_sample); %sort in ascending order
baseline_sorted_percentage=sortedVals(1:round(baseline_percentage*size(kymo_sample,1)),:); %take a specified percentage of these values
baseline_pc_median(k,:) = median(baseline_sorted_percentage); %find median of this percentage
end
baseline_pc_median = remv_zero_padding(baseline_pc_median, zeropad_len);
toc
end%% Gaussian derivative filter
function edge = Gaussian_Derivative_Filter_padded(data, thres)
tic;
disp('Applying 1st order Gaussian filter...')
% sigma for guassian filter
sigma_y = 2;
% length of zero-padding
zeropad_len = 50;
% zero-padding
edge = zero_padding(data, zeropad_len);
% 1st Order Derivative 1D Gaussian filter, detects slopes in time dimension
edge = Gaussian_Derivative_Filter(edge, sigma_y);
% remove zero-padding
edge = remv_zero_padding(edge, zeropad_len);
% correct for direction of the slope
edge = -edge;
% thresholding to remove slow flucturations in signal
edge(edge<thres)=0;
toc
end
function im = Gaussian_Derivative_Filter(im, sigma_y)
y_values = -ceil(4*sigma_y):ceil(4*sigma_y);
filter_y = calc1stOrderDerivative1DGaussian(y_values, sigma_y)';
im = imfilter(im, filter_y,'symmetric');
end
function [ D_values ] = calc1stOrderDerivative1DGaussian(x_values, sigma)
% Getting the 1D Gaussian values
G_values = calc1DGaussian(x_values, sigma);
% Obtaining the 1st order derivative of the
D_values = -2.*x_values./2./sigma^2.*G_values;
end
function [ G_values ] = calc1DGaussian(x_values, sigma)
G_values = 1/sigma/sqrt(2*pi).*exp(-x_values.^2./2./sigma^2);
end| Function | Description |
|---|---|
| imbinarize | Binarize 2-D grayscale image or 3-D volume by thresholding |
| graythresh | Global image threshold using Otsu's method |
| multithresh | Multilevel image thresholds using Otsu’s method |
| Function | Description |
|---|---|
| strel | Morphological structuring element |
| imdilate | Dilate image |
| imerode | Erode image |
| imopen | Perform morphological opening. Morphological opening is useful for removing small objects and thin lines from an image while preserving the shape and size of larger objects in the image. |
| imclose | Perform morphological closing. Morphological closing is useful for filling small holes in an image while preserving the shape and size of large holes and objects in the image. |
| bwskel | Skeletonize objects in a binary image. The process of skeletonization erodes all objects to centerlines without changing the essential structure of the objects, such as the existence of holes and branches. |
| bwperim | Find perimeter of objects in a binary image. A pixel is part of the perimeter if it is nonzero and it is connected to at least one zero-valued pixel. Therefore, edges of interior holes are considered part of the object perimeter. |
| bwhitmiss | Perform binary hit-miss transform. The hit-miss transform preserves pixels in a binary image whose neighborhoods match the shape of one structuring element and do not match the shape of a second disjoint structuring element.The hit-miss transforms can be used to detect patterns in an image. |
| imtophat | Perform a morphological top-hat transform. The top-hat transform opens an image, then subtracts the opened image from the original image.The top-hat transform can be used to enhance contrast in a grayscale image with nonuniform illumination. The transform can also isolate small bright objects in an image. |
| imbothat | Perform a morphological bottom-hat transform. The bottom-hat transform closes an image, then subtracts the original image from the closed image. The bottom-hat transform isolates pixels that are darker than other pixels in their neighborhood. Therefore, the transform can be used to find intensity troughs in a grayscale image. |
| bwmorph | Morphological operations on binary images. To perform morphological operations on a 3-D volumetric image, use bwmorph3. |
| Function | Description |
|---|---|
| imclearborder | Suppress light structures connected to image border |
| watershed | Watershed transform |
| Function | Description |
|---|---|
| regionprops3 | Measure properties of 3-D volumetric image regions. See also regionprops |
| bwconncomp | Find and count connected components in binary image |
| labelmatrix | Create label matrix from bwconncomp structure |
| bwlabeln | Label connected components in binary image |
| Function | Description |
|---|---|
| imtranslate | Translate image |
| imrotate3 | Rotate 3-D volumetric grayscale image |
| Function | Description |
|---|---|
| imregister | Intensity-based image registration |
| imregconfig | Configurations for intensity-based registration |
| imregdemons | Estimate displacement field that aligns two 2-D or 3-D images |
| imwarp | Apply geometric transformation to image |
| Function | Description |
|---|---|
| figure | Create figure window |
| subplot | Create axes in tiled positions |
| tiledlayout | Create tiled chart layout for displaying subplots |
| nexttile | Create axes in tiled chart layout |
| gca | Current axes or chart |
| gcf | Current figure handle |
| axis | Set axis limits and aspect ratios |
| xlim | Set or query x-axis limits. See also ylim and zlim |
| xlabel | Label x-axis. See also ylabel and zlabel |
| title | Add title |
| sgtitle | Add title to grid of plots |
| Function | Description |
|---|---|
| imagesc | Display image with scaled colors |
| contour | Contour plot of matrix. Good for showing boundary of ROIs in mask |
| plot | 2-D line plot |
| xline | Vertical line with constant x-value. See also yline |
| scatter | Scatter plot |
| imshowpair | Compare differences between images |
| stackedplot | Stacked plot of several variables with common x-axis |
| Function | Description |
|---|---|
| volshow | Display volume |
| orthosliceViewer | Browse orthogonal slices in grayscale or RGB volume |
| sliceViewer | Browse image slices |
| scatter3 | 3-D scatter plot |
| plot3 | 3-D point or line plot |
| view | Camera line of sight |
%% MAIN SCRIPT (or in command window)
% path to data
ops.filedir = '../originals/'; % folder
% path to saving directory
ops.savedir = '../outputs/'; % folder
if ~exist(ops.savedir, 'dir')
mkdir(ops.savedir)
end
loop_through_folder(ops);
disp('Done:)')
%--------------------------------%
% %
% (LOCAL) FUNCTIONS %
% %
%--------------------------------%
function loop_through_folder(ops)
%% loops through folder and subfolders
filelist = dir(foldername);
ops.savepath = ops.savedir;
for i = 1:length(filelist) % the first two are '.' and '..', skip
if strcmp(filelist(i).name,'.')
% % save data in the original data folder
% ops.savepath = filelist(i).folder;
% if ~exist(ops.savepath, 'dir')
% mkdir(ops.savepath)
% end
elseif strcmp(filelist(i).name,'..')
continue
elseif filelist(i).isdir % is a folder
disp(filelist(i).name);
% create new folder for saving data for each subfolder
ops.savedir = fullfile(ops.savedir, filelist(i).name);
if ~exist(ops.savedir, 'dir')
mkdir(ops.savedir)
end
ops.filedir = fullfile(filelist(i).folder,filelist(i).name);
loop_through_folder(ops);
elseif contains(filelist(i).name, '.czi') % change file format
ops.filename = fullfile(filelist(i).folder, filelist(i).name);
[~,filename,~] = fileparts(ops.filename);
% create new folder for saving data for each image file
ops.savedir = fullfile(ops.savepath, filename);
if ~exist(ops.savedir, 'dir')
mkdir(ops.savedir)
end
try
custom_function(ops);
catch
warning('Error. Check file.')
end
end
end
end
function your_processing(ops)
%--------------------------------%
% %
% CODES FOR PROCESSING %
% %
%--------------------------------%
tic;
[~,filename,~] = fileparts(ops.filename);
disp(filename);
% Put your tested pipeline here.
% You may remove the plots used for varifying the steps.
% It is recommended to save the key plots for debugging purpose.
% Add the line 'close(gcf)' to close the figure after saving.
toc;
endDownloading official Matlab add-ons to Clarke require admin right. Let Jane knows if you need any additional toolbox.
For other scripts that could be useful to your project, you may download them into /PATH/TO/PROJECT/readmes and add that to path using
addpath('/PATH/TO/PROJECT/readmes')For any github resources, it is a good idea to stay updated and keep track on the changes you've made to the scripts. Follow the step in here to set up set up source control.
Common lab SOPs:
Bioinformatics-related:
Image analyses-related:
Programming-related: