ImageJ / Fiji
Our team promotes the use of the open source platform ImageJ for bioimage processing and analysis.
With ImageJ, we develop custom solutions to process microscopy data, with mainly two development goals in mind: image analysis workflows to extract quantitative information from the images and tools to accelerate and automate image handling and data management.
Among the ImageJ macros we have developed, here are some example workflows addressing common problems. More recent macros and projects can be found in this respository: https://github.com/SebastienTs
Automated Multicellular Tissue Analysis
Segment the cells of a multicellular tissue in well contrasted images (Fig. 1.1). The geometry of the cells and their organization is automatically extracted and exported to an ImageJ results table. This includes: Cell area, major, minor fitted ellipse radii + major axis orientation, and number of neighbors. Manual correction of the segmented cells is optionally supported (merge / split objects). The analysis can also be restricted to a user defined region (Fig 1.2).
Fig. 1.1 Drosophila embryo expressing mCherry-cadherin (spinning disk, maximum intensity projection, courtesy of Annalisa Letiza, IBMB-CSIC)
Fig. 1.2 Segmentation results
1) Gaussian blur pre-filter radius: should be adapted to image noise, typical values in the range 1 to 3 pixels.
2) Cell detection sensitivity: values in range 0 to 255 (for 8-bit images). Can be set interactively (set to -1).
3) Minimum area overlap between dilated cell boundaries and neigbor cells. Used to discard neighbor cells only cornering a cell.
Measurements for each cell, the row index corresponds to the labels in Fig 1.2.
FISH Signals detection in Human Spermatozoids
Aim and prerequisites
Segments human spermatozoids nuclei from DAPI channel and classifies nuclei based on the number of FISH signals they contain. Reports the frequency occurence of user defined categories (combinations of spot multiplicity in FISH channels) as well as the position (point selections) of the detected nuclei falling in these categories.
The input image should be an hyperstack with 4 channels: DAPI (first channel) and three FISH channels. The images are typically obtained as a maximum intensity projection of a thin image stack.
The macro requires this LUT to run (should be copied to Fiji luts folder).
When launching the macro, the number of nucleus categories to detect can be set from the dialog box. The combination of FISH spots for each nucleus category is defined in the next dialog box. If a specific channel is un-ticked, the spot count is ignored for this channel.
Fig. 2.1 The original image is an hyperstack with 4 channels (DAPI first) and 3 FISH channels, represented here as a color composite.
(Image Courtesy of Anna Godo, Universitat Autònoma de Barcelona)
Nuclei segmentation parameters
1) Nucleus smooth: radius (pix) of the Laplacian of Gaussian filter used as a pre-filtering step.
2) Nucleus sensitivity: a negative value typically in the range [-0.5 to -0.1], the smaller the value the less sensitive the detection.
3) The next two values define the lower and upper bound for the nuclei area. These values are dimensionless multiplicative coefficients of the median area of all the particles detected in the image (isolated nuclei assumed in majority).
4) Nucleus minimum circularity from 0 (line) to 1 (perfect circle).
5) The number of dilation steps after nucleus segmentation before counting the spots (ensure spots close to nucleus edge are counted too).
Spots detection parameters
1) Spot smoothing: radius (pix) of the Laplacian of Gaussian pre-filter applied to each FISH channel.
2) Chan 1-3 spot level: detection sensitivity in each of the three FISH channels (8-bit scale).
The contours of the detected nuclei are color coded according to their spot content. Only the nuclei checking the user selected combination of spots are selected (point selection). Clusters of nuclei that could not be split are discarded. The frequency of the nuclei found in each of the user defined classes is reported.
Fig. 3.1 A fluorescence image from a protein microarray
(contrast saturated, Image courtesy of Fulvio Santacatterina, CBMSO)
1)The number of wells per row.
2) The number of wells per column.
3) Noise Tolerance: well detection sensitivity (the default value can be used if the original image is 16-bit with correct exposition).
4) Background subtraction can be used to compensate for background intensity bias.
The ROI labels correspond to the well labels overlaid on the image.
Blood vessel segmentation and network analysis
Aim and prerequisites
Segments blood vessels in a 3D stack. Suited for well contrasted images and reasonably uniform vessel width. The stack should be calibrated. This macro requires the ImageJ 3D suite that can be found here and this LUT (that should be copied to Fiji luts folder).
Fig. 4.1 3D rendering of a test image (lectin-rhodamine staining in mice, Macro SPIM acquisition)
(Image courtesy of Alexandre Calon, IRB Barcelona)
1) The first two parameters should both be set to the approximate expected radius of the blood vessels They respectively control tube closing (hollow vessels) and vessel pre-filter radius.
2) Vessel threshold: detection sensitivity (the lower the more sensitive).
3) Remove branches with end-points: skeleton branches leading to an end-point can optionally be trimed.
4) Dilate Skeleton for viewing: 3D viewer visualization only.
5) Number of threads: multi-threading, for best performance adjust this value to twice the number of cores of your CPU(s).
6) Save intermediate steps: Save all intermediate images to disk.
The results consist of:
1) A label mask of connected components (Fig 4.2).
2) A labelled skeleton with branching and end-points (Fig 4.3).
3) A log holding some statistics on the blood vessel network (Fig 4.4) such as the average blood vessel cross-section area and diameter, the total length of the network, the number of branches and the overall vascularization density (volume fraction occupied by the blood vessels).
Fig 4.2 Skeleton with end-points (white)
Fig 4.3 Connected components label mask and branching points (yellow)
Fig 4.4 Statistics displayed in Fiji log window
Trainable WEKA Segmentation batch processing
Aim and prerequisites
Batch processes all the images (TIFF and JPEG files) located in a user defined folder by calling Fiji Weka trainable segmentation to classify each pixel, and reports the areas of each class in a human readable results table. The classifier to be applied should be previously trained on a representative image and exported to file (Save classifier) to an empty folder named Results inside the folder with the images to be processed.
After launching the macro, pick any image from the folder to process.
Parameters (variables in the macro file)
WekaVersion: The macro has only be tested with Weka trainable segmentation v3.3.2 even though it should also work with a different version. The version should be updated in the first line of the macro, (string variable WekaVersion). To find out the version of Weka trainable segmentation that you run, call the plugin and copy the name of the plugin window to this string.
Classifier: The classifier is expected to have the default name (classifier.model), unless the string classifier is updated in the macro preamble.
1) A label mask per image with one gray level per class (Figure 9.2).
2) A results table (figure 9.3) reporting the image names and fractional class areas (one image per row).
Fig 9.1 Original image (RGB image, Trichrome)
(Image courtesy of Alexandre Calon, IRB Barcelona)
Figure 9.2 Label mask
background (black), collagen (orange), muscle (cyan), nuclei (white)
Fig 9.3 Results table reporting class fractional areas
Micro-tubule comets 2D tracking
Aim and prerequisites
Process spinning disk time-lapses to 1) segment the soma and reports its area (first frame only), then count the number of micro-tubule comets inside the soma across the time-lapse.
Manual ROI drawing (tick box): By default the macro segments the cell for all time points. If this fails it is possible to use a user drawn segmentation mask (then fixed across time). This mask is used to 1) restrict the comet detection and 2) report cell area.
Approximate radius: Comet approximate radius in pixels.
Detection threshold: Intensity sensitivity level for comet detection.
Border retraction: Cell boundary retraction (pix). Used to ignore bright edge signal (autofluorescence).
Min spot area: Minimum area of a comet (in pixel). Smaller comets particles are discarded.
Spot dilation: Dilation of the detected comets (in pixel). The dilation is used after comet detection to ensure an overlap between the same comet from frame to frame and avoid multiple counting of the same comet.
1) A reference image showing the detected soma after border retraction and comets (Fig 10.2).
2) A log showing the image name, the soma area and the comet (spot) count.
Fig 10.1 A time point of the original time-lapse
Fig 10.2 The same time point with segmented soma (after border retraction) and segmented comets
Aim / Results
Attenuate close to parallel stripes, such as the ones typically observed in light sheet microscopy. The filter is designed to cope with some stripe angular spread (e.g. due to light sheet refraction at the sample surface as in Fig 5.1), and it can process a 3D stack (processing performed slice by slice).
Alternatively check this other project: https://github.com/SebastienTs/FFTdeStriper
Fig. 5.1 Original image (Zebrafish epiboly)
(Image courtesy of Maria Marsal, IBMB - CSIC)
Fig. 5.2 Filtered image with default settings
For stripe width similar to the sample image, setting "SubRad" and "OpenRad" to 40 pix and disabling sharpening (SharpenRad = 0) is a good starting point. Scale these values proportionally to the typical width of the stripes. Some guidelines to further adjust the parameters:
1) Thin stripes apparent in filtered image --> decrease OpenRad
2) Large, uniform, stripes apparent in filtered image --> increase SubRad
3) Filtered image looks "blocky" and too smooth --> increase OpenRad
4) Filtered image looks dirty and distorted --> decrease SubRad
SharpenRad / SharpenFeedback are the parameters of a post-processing high-pass filter meant to boost the high frequency part of the spectrum that might have been too attenuated by the main filter. Increase SharpenFeedback to sharpen more (practical range: 0 to 0.5), adjust SharpenRad to obtain the best results.
Unfold a tubular structure and flatten its surface (like peeling of and flattening the skin of a banana). Technically, compute the radial average intensity projection inside a ring centered on the radial symmetry axis of the object. The final image is a radial mapping of the intensity (radial angle along X, axial length along Y). Process only a single channel 3D stack but it is easy to process multiple channels by exporting and importing ROI manager selections (so as to perform the same operation).
The tube axis and contour are user drawn from three different views. The contours are saved to the ROI manager automatically. If the macro is launched with exactly 3 selections in the ROI manager it will use these contours automatically without asking the user to draw them (same transformation applied); it is hence easy to process multi-channel 3D stacks. The contours can also be saved from the ROI manager window to keep them as reference with the original data.
Fig. 6.1 The original 3D stack (composite 3D rendering of the 3 channels)
(Image Courtesy of Nareg Djabrayan, IRB Barcelona)
Fig. 6.2 Resulting unfolded image. The image has also been circularly shifted to display the best orientation
Somata segmentation and time response analysis in acute rodent brain slices
Detect cell somata in acute rodent brain slices time-lapses after noise reduction and image registration (as required), extract fluorescence changes relative to baseline (dF/F) over time, and classify cells as responsive versus non-responsive based on the presence or absence of significant fluorescence peaks in their time response. ImageJ version >=1.52 required.
Fig. 9.1 (1) Original time-lapse (filtered + maximum intensity projection), (2) somata detection, (3) normalized time traces extraction showing a responsive cell (top) and a non-responsive cell (bottom), (4) cell classification (responsive: red, non-responsive: yellow)
(Image courtesy of Hyojung Lee, IBEC)
Steps implemented and parameters (set in the macro code, not from a dialog box)
I) Noise reduction of the OGB-1 fluorescence intensity image sequences by:
1. Spatial filtering (Gaussian blur sigma: FltRad)
2. Temporal averaging (optional; TmpGroup: 1 or 2 frames averaging were used for this study)
II) Registration (image stabilization) if the sample suffers from spatial drift during the acquisition (set Register to true).
III) Identification of cell somata as intensity regional maxima in the maximum intensity temporal projection of the filtered time-lapse. Users can finely adjust somata detection by:
1. Adjusting noise tolerance (NoiseTol)
2. Interactively editing the detected somata from the ROI Manager
IV) Measurement of average intensity as a function of time inside fixed-radius circular regions (adjustable radius CellRad; typically 5-9 pixels for this study,) centered on detected maxima.
V) Independent normalization of time traces by scaling them by their average over a selected number of initial time frames (BasalRange frames of the spatially filtered movie included).
VI) Classification of the somata as responsive (red) or non-responsive (yellow) based on the presence/absence of significant peaks (typically 20%-30% increase respect to baseline, set by Alpha) in their normalized time trace. A cell is classified as positive if at least MinPeaks of minimum width MinWidth are detected.
VII) Localization of burst maxima / width (default at half height, adjustable by PeakEdthFraction) in the time traces.
VIII) Exportation of all results as tabular files: raw time traces, estimated initial basal levels, and detected bursts.
Note: For debugging purpose it is possible to sequentially plot all time traces and display detected peaks by setting PlotEach to true.
(Z,T) Hyperstack Stitcher
Aim / Results
Builds a stitched image from a muti-position 3D + time hyperstack. The XY positions of the montage should be coded as channels in the input hyperstack. Channel ordering can be configured in the dialog box to adapt to Column/Row and Meander/Comb configurations: The images should appear in this order when browsing the hyperstack with the channel slider. Fine stitching is supported (requires sufficient overlap between the views). The XY displacements of each field of view for stitching are computed for a single reference (Z,T) slice (user configurable) and applied to all slices (Z and T).
Note: for OME file naming the hyperstack can easily be generated by importing the files and transforming the stack to an hyperstack with xyzct convention (provided the L and T fields of the OME files are equal) or xyztc convention if the L field is constant. In this case the Column/column comb configuration should be used regardless of the configuration of the scan acquisition. In case of different file naming (e.g. single index) another convention might be required.
Fig. 7.1 Example stitching for 2x2 XY positions
(Image courtesy: Jelena Urosevic, IRB)
1) The number of Y positions of the acquisition (number of rows).
2) The number of X positions of the acquisition (number of columns).
3) The order in which the images appear when browsed from the hyperstack (channel slider).
4) The overlap between the images in both X and Y dimensions (coarse estimate).
5) A parameter used to control the tolerance of the stitiching.
6) Tick "stitch" only if you want to perform a fine stitching of the images, if not the coarse overlap estimate is used to position the images.
7) The reference z slice to perform the fine stitching (same X/Y displacements are applied to all other slice).
Aim / Results
Stitch a (Z,T,C) data set. No limit is set on the number of Z slices and time frames. The input is a folder with the raw tiff images (one image per file) as typically exported by motorized microscopes. These files must all be stores in the same folder and the file naming should ideally comply to OME-TIFF. Only --X, --Y and --Z fields with user defined number of digits are compulsory. --T, --C and --L fields with user defined number of digits are necessary for multiple time frames / channels data sets. A compatible data set is provided as a .zip archive. Before processing it unzip it to a given location.
The stitching is performed in a reference Z slice (and in a specific reference time frame and channel). The same offsets are applied to all the Z slices, time frames and channels. Before starting the batch processing a montage with the original images of the selected Z slice / time frame / channel is displayed together with the stitched image in this stack. If you are not satisfied with the result you can select another reference. The stitching is then performed time frame by time frame and slice by slice and the stitched images are exported to a single user defined output folder.
The code can also process a data set with multiple channels, the stitching is then computed once on a reference channel and then applied to the other channels.
1) The number of fields of view per column to stitch (not necessarily all).
2) The number of fields of view per row to stitch (not necessarily all).
3) The number of Z slice to process (not necessarily all).
4) The number of time frames that you want to process (not all).
5) The approximate overlap (%) between the images in both X and Y dimensions.
6) A parameter used to control the tolerance of the stitching: 0 (loose) to 1 (picky)
7) The image fusion mode you want to apply for the stitching
8) The common root name of the stitched image exported to the result folder
1) You can display the accumulated Z maximum intensity projection as progressing through the stack.
2) In case you want to process another channel tick "Additional channel". You then have to input the L and C fields indexes of the new channel. The processing can be performed by launching the macro N times (for N channels). In case you want to batch process two channels in a row you have to tick "Launch additional channel on complete", select one of the two channels as reference and configure the other channel in the C/L fields.
Aim and prerequisites
Experimental. Enable to interact with a large, single channel, z-stack (possibly exceeding the main memory of the computer) and to extract a volume of interest by marking several reference points. The code requires this plugin that should be copied to ImageJ plugins folder. It is only compatible with 8-bit stacks. The 3D stack should be opened first. The maximum intensity projection of the stack is computed from multiple angles. If the 3D stack is very large and does not fit into memory in should be opened as virtual stack.
Parameters and mode of operation
1) Step sets the angular step in degrees between the views computed in the maximum intensity projections stack.
2) Max angle sets the maximum (last) angle to compute.
3) Scale factor is a downscaling factor. Should be set so that the downscaled stack is less than half the main memory available. The estimated size of the downscaled stack is printed to screen, you can abort the macro by pressing escape if this size is too large.
4) Anaglyph can be ticked to generate a Red/Cyan anaglyph
5) In case the anaglyph is generated the angle step multiplier sets the offset between the views. It is always a multiple of the projection stack angular step.
Interactively add reference points to the projection stack by holding "shift" and left clicking. The 3D position of the reference points are reported to the log window (pixel unit in the original stack). If the points do not lock in well, use "alt" and left click instead when the feature of interest is "in front"..
Press "space" after adding reference points. At least 3 reference points should be set, the smallest volume of interest containing them is extracted (a new stack is created) and displayed in the 3D viewer.
Note: To extract several VOI from the same stack the maximum intensity projection stack and the depth projection stack can be saved to file and the macro launched again with the 3 stacks opened (original + projection stacks). The projections will then not be recomputed.
Fig. 8.1 The downscaled interactive intensity projection stack
(Image courtesy of Alexandre Calon, IRB Barcelona)
Fig. 8.2 The extracted Volume of Interest
Track loosely convex moving objects (possibly dividing) in a binary stack and label them with a unique ID. Optionally output a division file compatible with Cell Tracking Challenge.
ErodRad: Erosion radius (increase to avoid false linking but objects with smaller radius are lost).
NoiseTol: Noise tolerance to split out touching objects (0: sensitive, higher: less sensitive).
DivFile: Path to the division text fil, leave empty not to export.