ImageJ

{{scripting}}
== Introduction ==

This page is a primer of '''ImageJ2 only''' Python scripts. It means that the examples included here avoid IJ1 as much as possible, unless it's really necessary.

== Scripts ==

=== Z Projection ===

<source lang='python'>
# @ImageJ ij
# @Dataset data
# @String(label="Projection Type",choices={"Max","Mean","Median","Min", "StdDev", "Sum"}) proj_type


from net.imagej.axis import Axes
from net.imagej.ops import Ops


def proj_method(method):
    return {
        'Max': Ops.Stats.Max,
        'Mean': Ops.Stats.Mean,
        'Median': Ops.Stats.Median,
        'Min': Ops.Stats.Min,
        'StdDev': Ops.Stats.StdDev,
        'Sum': Ops.Stats.Sum,
    }.get(method, Ops.Stats.Max)


def main():

    # Select which dimension to project
    z_dim = data.dimensionIndex(Axes.Z)

    if z_dim == -1:
        print("Z dimension not found.")
        return

    if data.dimension(z_dim) < 2:
        print("Z dimension has only one frame.")
        return

    # Write the output dimensions
    projected_dimensions = [data.dimension(d) for d in range(0, data.numDimensions()) if d != z_dim]

    # Create the output image
    z_projected = ij.op().create().img(projected_dimensions)

    # Create the op and run it
    proj_op = ij.op().op(proj_method(proj_type), data)
    ij.op().transform().project(z_projected, data, proj_op, z_dim)

    # Create a dataset
    z_projected = ij.dataset().create(z_projected)

    # Set the correct axes (is that needed ?)
    axes = [data.axis(d) for d in range(0, data.numDimensions()) if d != z_dim]
    z_projected.setAxes(axes)

    print(z_projected)
    ij.ui().show("z_projected", z_projected)

main()

</source>

=== Apply Threshold ===


<source lang='python'>
# @ImageJ ij
# @Dataset data

method_threshold = "otsu"

# Get the histogram
histo = ij.op().run("image.histogram", data.getImgPlus())

# Get the threshold
threshold = ij.op().run("threshold.%s" % method_threshold, histo)

# Apply the threshold
thresholded = ij.op().run("threshold.apply", data.getImgPlus(), threshold)

# Create output
thresholded = ij.dataset().create(thresholded)

ij.ui().show("thresholded", thresholded))
</source>

A more direct way if you don't need to modify the threshold is:

<source lang='python'>

# @ImageJ ij
# @Dataset data

method_threshold = "otsu"

# Apply the threshold
thresholded = ij.op().run("threshold.%s" % method_threshold, data.getImgPlus())

# Create output
thresholded = ij.dataset().create(thresholded)

ij.ui().show("thresholded", thresholded)
</source>


=== Subtract stack to its first image ===

<source lang='python'>

# @ImageJ ij
# @Dataset data

from net.imglib2.util import Intervals
from net.imagej.axis import Axes
from net.imagej.ops import Ops

# Convert input
converted = ij.op().convert().float32(data.getImgPlus())

# Get the first frame (TODO: find a faser way !)
t_dim = data.dimensionIndex(Axes.TIME)
interval_start = []
interval_end = []
for d in range(0, data.numDimensions()):
	if d != t_dim:
		interval_start.append(0)
		interval_end.append(data.dimension(d) - 1)
	else:
		interval_start.append(0)
		interval_end.append(0)
		
intervals = interval_start + interval_end
intervals = Intervals.createMinMax(*intervals)

first_frame = ij.op().transform().crop(converted, intervals)

# Allocate output memory (wait for hybrid CF version of slice)
subtracted = ij.op().create().img(converted)

# Create the op
sub_op =  ij.op().op("math.subtract", first_frame, first_frame)

# Setup the fixed axis
fixed_axis = [d for d in range(0, data.numDimensions()) if d != t_dim]

# Run the op
ij.op().slice(subtracted, converted, sub_op, fixed_axis)

# Clip image to the input type
clipped = ij.op().create().img(subtracted, data.getImgPlus().firstElement())
clip_op = ij.op().op("convert.clip", data.getImgPlus().firstElement(), subtracted.firstElement())
ij.op().convert().imageType(clipped, subtracted, clip_op)

# Show it
ij.ui().show("subtracted", clipped)
</source>


=== Apply DOG Filter ===

<source lang='python'>
# @ImageJ ij
# @Dataset data

from net.imagej.axis import Axes

sigma1 = 4.2
sigma2 = 1.25

# IJ Ops version

pixel_type = data.getImgPlus().firstElement().class
converted = ij.op().convert().float32(data.getImgPlus())

# Allocate output memory (wait for hybrid CF version of slice)
dog = ij.op().create().img(converted)

# Create the op
dog_op = ij.op().op("filter.dog", converted, sigma1, sigma2)

# Setup the fixed axis
t_dim = data.dimensionIndex(Axes.TIME)
fixed_axis = [d for d in range(0, data.numDimensions()) if d != t_dim]

# Run the op
ij.op().slice(dog, converted, dog_op, fixed_axis)

# Clip image to the input type
clipped = ij.op().create().img(dog, data.getImgPlus().firstElement())
clip_op = ij.op().op("convert.clip", data.getImgPlus().firstElement(), dog.firstElement())
ij.op().convert().imageType(clipped, dog, clip_op)

# Show it
ij.ui().show("dog", clipped)
</source>

=== Apply a mask ===

Given a mask (binary image) and a raw image, remove background pixel from raw by keeping only those in the mase (different from 0).

<source lang='python'>

# @ImageJ ij
# @Dataset data
# @Dataset mask

output = ij.dataset().create(data)

targetCursor = output.localizingCursor()
dataRA = data.randomAccess()
maskRA = mask.randomAccess()

while targetCursor.hasNext():
	targetCursor.fwd()
	dataRA.setPosition(targetCursor)
	maskRA.setPosition(targetCursor)

	if maskRA.get().get() == 0:
		targetCursor.get().set(0)
	else:
		targetCursor.get().set(dataRA.get())

ij.ui().show("output", output)
</source>

As specified by @stelfrich on Gitter, the particular case when foreground pixel are 1 and background pixels are 0 can be simpler to write with a multiplication of the two images.

<source lang='python'>
# @ImageJ ij
# @Dataset data
# @Dataset mask
 
output = ij.dataset().create(data)
ij.op().op("math.multiply", output, data.getImgPlus(), mask.getImgPlus())
ij.ui().show("output", output)
</source>

=== Retrieve objects/particles from a mask ===

Get a list of <code>LabelRegion</code> and display the center of mass with the IJ1 ROI Manager.

<source lang='python'>
# @ImageJ ij
# @Dataset data
# @Dataset mask

from net.imagej.axis import Axes

from net.imglib2.algorithm.labeling.ConnectedComponents import StructuringElement
from net.imglib2.roi.labeling import LabelRegions

from ij.gui import PointRoi
from ij.plugin.frame import RoiManager

def get_roi_manager(new=False):
	rm = RoiManager.getInstance()
	if not rm:
		rm = RoiManager()
	if new:
		rm.runCommand("Reset")
	return rm

img = mask.getImgPlus()

labeled_img = ij.op().run("cca", img, StructuringElement.EIGHT_CONNECTED)

regions = LabelRegions(labeled_img)
region_labels = list(regions.getExistingLabels())

print("%i regions/particles detected" % len(region_labels))

# Now use IJ1 RoiManager to display the detected regions

rm = get_roi_manager(new=True)
for label in region_labels:
 
    region = regions.getLabelRegion(label)
 
    center = region.getCenterOfMass()
    x = center.getDoublePosition(0)
    y = center.getDoublePosition(1)
 
    roi = PointRoi(x, y)
    if center.numDimensions() >= 3:
        z = center.getDoublePosition(2)
        roi.setPosition(int(z))
     
    rm.addRoi(roi)

    # You can also iterate over the `data` pixel by LabelRegion

	cursor = region.localizingCursor()
	dataRA = data.randomAccess()
	while cursor.hasNext():
		cursor.fwd()
		dataRA.setPosition(cursor)
	
		x = cursor.getDoublePosition(0)
		y = cursor.getDoublePosition(1)

		# Pixel of `data`
		pixel = dataRA.get()

		# Do whatever you want here
		# print(x, y, pixel)
</source>

=== Extract patch ===

NOT OPTIMIZED - WIP

Extract a patch from an image given its center, its size and also an orientation.

<source lang='python'>
# @ImageJ ij
# @Dataset data

from net.imglib2.interpolation.randomaccess import LanczosInterpolatorFactory
from net.imglib2.view import Views
from net.imagej.axis import Axes

import math


def extract_patch(data, origin, orientation, size):
	"""Inspired from A. Matov and al., 2010, Nature Methods
	"""

	e1 = [-math.cos((orientation + 90) * math.pi / 180), math.sin((orientation + 90) * math.pi / 180), 0]
	e2 = [math.sin((orientation + 90) * math.pi / 180), math.cos((orientation + 90) * math.pi / 180), 0]
	e3 = [0, 0, 1]
	v1 = [-size//2, size//2]
	v2 = [-size//2, size//2]
	v3 = [0, 0]

	# Generate a list of coordinates
	listOfCoords = []
	for z in range(v3[0], v3[1] + 1):
		for y in range(v2[0], v2[1] + 1):
			for x in range(v1[0], v1[1] + 1):
				coord = []
				coord.append(x * e1[0] + y * e2[0] + z * e3[0] + origin[0])
				coord.append(x * e1[1] + y * e2[1] + z * e3[1] + origin[1])
				coord.append(x * e1[2] + y * e2[2] + z * e3[2] + origin[2])
				listOfCoords.append(coord)
	
	
	# Create realRandomAccess cursor
	interpolator = LanczosInterpolatorFactory()
	interpolant = Views.interpolate(Views.extendZero(data.getImgPlus()), interpolator)
	realRandomAccess = interpolant.realRandomAccess()
	
	# Create the output
	patch_dims = [v1[1] - v1[0] + 1, v2[1] - v2[0] + 1, v3[1] - v3[0] + 1]
	# TODO : setup same Axes as input
	axes = [Axes.X, Axes.Y, Axes.Z]
	patch = ij.dataset().create(patch_dims, "patch", axes, data.getValidBits(), data.isSigned(), False)
	
	patchRA = patch.randomAccess()
	
	i = 0
	for x in range(patch_dims[0]):
		for y in range(patch_dims[1]):
			for z in range(patch_dims[2]):
			
				patchRA.setPosition([x, y, z])
				realRandomAccess.setPosition(listOfCoords[i])		
				patchRA.get().set(realRandomAccess.get())
				i += 1
				
	return patch

# Extract a patch of size 20 pixel at a defined position with a defined orientation
orientation = 46
origin = [294, 81, 0]
size = 10
patch = extract_patch(data, origin, orientation, size)
ij.ui().show(patch)
</source>

=== Crop an image along one single axis ===

<source lang='python'>
from net.imagej.axis import Axes
from net.imglib2.util import Intervals


def get_axis(axis_type):
    return {
        'X': Axes.X,
        'Y': Axes.Y,
        'Z': Axes.Z,
        'TIME': Axes.TIME,
        'CHANNEL': Axes.CHANNEL,
    }.get(axis_type, Axes.Z)


def crop_along_one_axis(ij, data, intervals, axis_type):
	"""Crop along a single axis using Views.

	Parameters
	----------
	intervals : List with two values specifying the start and the end of the interval.
	axis_type : Along which axis to crop. Can be ["X", "Y", "Z", "TIME", "CHANNEL"]
	"""

	axis = get_axis(axis_type)
	interval_start = [data.min(d) if d != data.dimensionIndex(axis) else intervals[0] for d in range(0, data.numDimensions())]
	interval_end = [data.max(d) if d != data.dimensionIndex(axis) else intervals[1] for d in range(0, data.numDimensions())]

	interval = interval_start + interval_end
	interval = Intervals.createMinMax(*interval)

	output = ij.op().run("transform.crop", data, interval, True)

	return output
</source>

=== Rotate all the frames of a stack ===

<source lang='python'>
# @ImageJ ij
# @Dataset data

import math

from net.imagej.axis import Axes
from net.imglib2.util import Intervals

from net.imglib2.realtransform import RealViews
from net.imglib2.realtransform import AffineTransform2D
from net.imglib2.interpolation.randomaccess import LanczosInterpolatorFactory
from net.imglib2.view import Views


def get_axis(axis_type):
    return {
        'X': Axes.X,
        'Y': Axes.Y,
        'Z': Axes.Z,
        'TIME': Axes.TIME,
        'CHANNEL': Axes.CHANNEL,
    }.get(axis_type, Axes.Z)
    

def crop_along_one_axis(ij, data, intervals, axis_type):
	"""Crop along a single axis using Views.

	Parameters
	----------
	intervals : List with two values specifying the start and the end of the interval.
	axis_type : Along which axis to crop. Can be ["X", "Y", "Z", "TIME", "CHANNEL"]
	"""

	axis = get_axis(axis_type)
	interval_start = [data.min(d) if d != data.dimensionIndex(axis) else intervals[0] for d in range(0, data.numDimensions())]
	interval_end = [data.max(d) if d != data.dimensionIndex(axis) else intervals[1] for d in range(0, data.numDimensions())]

	interval = interval_start + interval_end
	interval = Intervals.createMinMax(*interval)

	output = ij.op().run("transform.crop", data, interval, True)

	return output

angle = 90

# Get the center of the images so we do the rotation according to it
center = [int(round((data.max(d) / 2 + 1))) for d in range(0, data.numDimensions())]

# Convert angles to radians
angle_rad = angle * math.pi / 180

# Build the affine transformation
affine = AffineTransform2D()
affine.translate([-p for p in center])
affine.rotate(angle_rad)
affine.translate(center)

# Get the interpolator
interpolator = LanczosInterpolatorFactory()

# Iterate over all frame in the stack
axis = Axes.TIME
output = []
for d in range(data.dimension(axis)):

	# Get the current frame
	frame = crop_along_one_axis(ij, data, [d, d], "TIME")

	# Get the interpolate view of the frame
	extended = ij.op().run("transform.extendZeroView", frame)
	interpolant = ij.op().run("transform.interpolateView", extended, interpolator)

	# Apply the transformation to it
	rotated = RealViews.affine(interpolant, affine)
	
	# Set the intervals
	rotated = ij.op().transform().offset(rotated, frame)

	output.append(rotated)

output = Views.stack(output)

ij.ui().show("output", output)
</source>

== Resources ==
* [[Jython_Scripting|Jython Scripting]]
* [[Jython_Scripting_Examples|Jython Examples]]

[[Category:Scripting]]
[[Category:ImageJ2]]
[[Category:Ops]]
[[Category:SciJava]]