This notebook dives into some technical details of the ImageJ Ops library.
It is recommended that you first read and understand the ImageJ Ops notebook before tackling this one.
%classpath config resolver scijava.public "https://maven.scijava.org/content/groups/public"
%classpath add mvn net.imagej imagej 2.0.0-rc-71
ij = new net.imagej.ImageJ()
ij.getVersion()
Added new repo: scijava.public
2.0.0-rc-71
This tutorial covers the so-called SPECIAL
ops. Every special op produces one primary output from some fixed number of primary inputs. The number of such inputs for a given op is known as its ARITY
. The framework provides interfaces for three arities:
Additional arities (e.g., ternary) are feasible, but not implemented.
The SPECIAL
ops can be classified into three main categories: COMPUTER
, FUNCTION
and INPLACE
:
COMPUTER
: An op which will write its result to pre-allocated output, i.e. an empty (zero) imageFUNCTION
: An op which will allocate on its own and then return itINPLACE
: An op which will write its output to one of the input arguments, hence mutating one of the inputsSpecial ops may also have any number of additional secondary inputs, which do not affect the arity of the op. More on that below.
Let's try some example op calls, to illustrate the difference between these various kinds of special ops.
First, we will create some test images to work with. In order to do that, we define and use a helper function, which creates an image of specified width and height using the create.img
op, and then populates the contents of that image using the image.equation
op, which sets the value of each sample according to a specified formula. In this case, we create four images depicting different patterns based on the location of a sample, which may be accessed via the variable p[]
inside the formula.
import net.imglib2.type.numeric.real.DoubleType
// Helper function to create DoubleType images from a formula
imageFromFormula = { name, width, height, formula ->
img = ij.op().run("create.img", [width, height], new DoubleType())
ij.op().run("image.equation", img, formula)
println("Created image: " + name + " of size [" + width + "," + height + "]")
return img
}
sinusoid = imageFromFormula("Sinusoid", 150, 150, "10 * (Math.cos(0.3*p[0]) + Math.sin(0.5*p[1]))")
gradient = imageFromFormula("Gradient", 150, 150, "p[0] + p[1]")
spot = imageFromFormula("Circle", 150, 150, "Math.sqrt((p[0]-75)*(p[0]-75) + (p[1]-75)*(p[1]-75))")
cross = imageFromFormula("Diamond", 150, 150, "Math.abs(p[0]-75) + Math.abs(p[1]-75)")
ij.notebook().display([["sinusoid": sinusoid, "gradient": gradient, "spot": spot, "cross": cross]])
Created image: Sinusoid of size [150,150] Created image: Gradient of size [150,150] Created image: Circle of size [150,150] Created image: Diamond of size [150,150]
sinusoid | gradient | spot | cross |
---|---|---|---|
First we will look at binary ops. A BINARY
op is one with two primary inputs, which produces an output. Suppose we want to add two images together. There are several versions of the math.add
op. In the following example we will try the functional version:
//Pointwise add FUNCTION
functionOutput = ij.op().run("math.add", sinusoid, gradient)
ij.notebook().display(functionOutput)
The above invocation calls a math.add
FUNCTION
, which allocates and returns a newly created output image, without mutating any of the input data. This is great from a functional standpoint, but costs space and time. Now suppose we already have a properly dimensioned output image allocated and ready to go. In that case, it would be wasteful to call math.add as a function, only to then copy the data from the newly allocated result back into our already-existing image. Much better would be to tell math.add
to use our existing output buffer directly. Fortunately, this approach is made possible by COMPUTER
ops, which require you to specify a pre-allocated output buffer where the result will be stored.
For example, imagine you want to do some analysis on a large set of images which would require to run e.g. a Gaussian smoothing as an intermediate operation. In this case, this intermediate step would not have to be store as it is not part of the end result. Thus, we can make use of the COMPUTER
paradigm. We pre-allocate a single output buffer to store the results of the smoothing and continue with our analysis for the current image. For the next image, we simply overwrite the same buffer and continue. This approach is much faster compared to calculating the smoothing using a FUNCTION
which would create a new output buffer for every intermediate step.
In the following example, we will first pre-allocate an output image and then write the result of the addition to this image:
// Perallocate output image
computerOutput = ij.op().run("create.img", sinusoid)
// Pointwise add using a COMPUTER
ij.op()run("math.add", computerOutput, spot, gradient)
ij.notebook().display(computerOutput)
In the above call, the same operation is performed, but the output is written into the provided "result" buffer. COMPUTER
ops do have an important restriction: the output buffer MUST NOT be the same reference as any of the inputs. Otherwise, the computer would end up mutating the input data as a consequence of writing to the output buffer. Hence, computers cannot be used to calculate results "in place" -- that's what INPLACE
ops are for.
import net.imagej.ops.special.inplace.Inplaces;
import net.imagej.ops.Ops
// Create a copy of one of the example images which we want to mutate
gradientCopy = ij.op().run("copy.rai", gradient)
// Create an INPLACE math.add Op
inplaceOp = Inplaces.binary1(ij.op(), Ops.Math.Add.class, gradientCopy, sinusoid)
""
inplaceOp.mutate(gradientCopy)
max = ij.op()run("stats.max", gradientCopy)
println("The maximum pixel value is: " + max)
ij.notebook().display(gradientCopy)
The maximum pixel value is: 308.6507787063433
The above invocation differs from earlier: the Inplaces.binary1
method requests a math.add
binary op which mutates its first argument inplace. A BinaryInplace1Op instance is returned, upon which we then call mutate(). In order to show that the content of the mutated image actually changes, we also calculate and print the maximum pixel value of the image using the stats.max
op. If the above cell is executed more than once, one can see that the maximum of the image actually changes as values are added to the pixels. Note that we stored the INPLACE
op in a variable. Such ops can be used in other ways too; see REUSING SPECIAL OPS
below.
In addition to binary ops which have the ARITY
2, there are also UNARY
(one primary input, one primary output) and NULLARY
(no primary inputs, one primary output). Next, we are going to look at an example for both. First, we look at the UNARY
math.sqrt
op which simply calculates the square root of a number. Note, that this op is also a COMPUTER
, meaning we first create an output container to write the result to.
import net.imglib2.type.numeric.real.DoubleType
// Create preallocated output
number = new DoubleType();
// Calculate the square root of 64
ij.op().run("math.sqrt", number, new DoubleType(64))
"The sqrt of " + 64 + " is " + number
The sqrt of 64 is 8.0
Second, we look at the NULLARY
math.zero
op. This op simply sets a number to zero. Note, that the variable number
still contains the square root of the calculation from the cell above.
import net.imglib2.type.numeric.real.DoubleType
ij.op().run("math.zero", number)
"number == 0.0? : " + number.equals(new DoubleType(0))
number == 0.0? : true
One of the most useful characteristics of special ops is the ability to look them up in a type-safe way, obtaining a reference to an op instance which can then be used however you need -- one time or many.
import net.imagej.ops.special.computer.Computers
import net.imagej.ops.Ops
// Create output for the COMPUTER
computerOutput = ij.op().run("create.img", sinusoid)
// Get an instance of the "math.add" op and store it for later
addOp = Computers.binary(ij.op(), Ops.Math.Add.class, computerOutput, gradient, cross)
net.imagej.ops.math.IIToIIOutputII$Add@3ce19b95
In this case, we use the Computers.binary
utility method to require a BINARY
COMPUTER
Op
, so that we can take advantage of the BinaryComputerOp API, which is type-safe -- whereas the Op interface's run("op.name", inputs..)
method is not.
As usual, compile-time type safety is a two-edged sword: helpful to avoid coding errors, but sometimes rather verbose. As this notebook uses Groovy, which is a Scripting language, we can get around some of Java's verboseness. In java, the line from above creating the op would look like this:
final BinaryComputerOp<IterableInterval<DoubleType>, IterableInterval<DoubleType>, IterableInterval<DoubleType>> addOp = Computers.binary(ops, Ops.Math.Add.class, computerOutput, gradient, cross);
However, the verbose way of doing it helps to understand what it going on underneath.
In the cell below we can now just call the run()
method to execute the op using the inputs specified above.
addOp.run()
ij.notebook().display(computerOutput)
Or we just run the op specifying different arguments by calling the compute()
method. Note that for the compute()
method the order of arguments will change. In this case the output to write to will be the last argument. Also note that the name of the method coincides with the type of op we are using. Here, we are using a COMPUTER
op, hence the method is called compute()
. In the case of an INPLACE
op, the method will be called mutate()
as one of the inputs is mutated inplace.
// Call the op using different inputs
addOp.compute(gradient, spot, computerOutput)
ij.notebook().display(computerOutput)
An op may additionally have secondary inputs. These are usually algorithm parameters or an ENUM
that specifies the behaviour of the op. One example would be the sigma value or the out-of-bounds strategy used for a Gaussian filter op, which would look like this (sigma = 2, out-of-bounds = Boundary.SINGLE):
gaussSigma2ComputerOp = Computers.unary(ij.op(), Ops.Filter.Gauss.class, input, input, 2, new OutOfBoundsMirrorFactory(Boundary.SINGLE))
Note the difference between primary and secondary inputs: the latter are not given when calling the special op's type-safe API; hence, secondary inputs remain fixed as given during the original matching (when we requested the op). Thus, the UNARY
COMPUTER
we saved in the variable gaussSigma2ComputerOp
will always perform a Gaussian filtering with sigma 2. Example usage could look like this:
gaussSigma2ComputerOp.compute(image, smoothedImage)
The variable smoothedImage
will now contain the smoothed image using a sigma value of 2 as defined above when we requested the op.
Which inputs are considered "primary" is a decision left to the author of each op. But here are two rules of thumb:
It is important to understand that each call into the Ops framework requires the system to search for the best-matching op amongst all those available. This is very useful for reasons of extensibility: the best op for the job is always selected, and anyone can extend the system with new ops optimized for particular scenarios. However, matching can have substantial performance implications:
import net.imglib2.type.numeric.real.DoubleType
import net.imglib2.IterableInterval
slowStart = System.currentTimeMillis();
200.times {
// Search for the op every time and then run it
ij.op().run("math.add", computerOutput, gradient, spot)
}
slowEnd = System.currentTimeMillis();
println("ij.op().run() needs\t" + (slowEnd - slowStart) + " ms for 200 iterations")
fastStart = System.currentTimeMillis();
200.times {
// Execute the same op instance repeatedly
addOp.compute(gradient, spot, computerOutput)
}
fastEnd = System.currentTimeMillis();
println("addOp.compute() needs \t" + (fastEnd - fastStart) + " ms for 200 iterations")
""
ij.op().run() needs 2812 ms for 200 iterations addOp.compute() needs 14 ms for 200 iterations
As one can see, reusing the same matched op every time is much faster. However, it is the responsibility of the caller to ensure that subsequent arguments passed in are compatible with that same op.
Some ops accept other ops as arguments, forming chains or trees of ops. For example, the map op executes the given op over all elements of an iteration. In the following example, we want to take the pixelwise square root of an image using the math.sqrt
op. Unfortunately, this op is only applicable for singe numbers compared to the math.add
op which can already handle IterableInterval
for convenience. However, we can make use of the map
op which expects another op as input and executes it on every element of an iterable; in this case every pixel.
import net.imagej.ops.special.computer.Computers
import net.imglib2.type.numeric.real.DoubleType
import net.imagej.ops.Ops
// Create output for the COMPUTER
mapOutput = ij.op().run("create.img", sinusoid)
imgMax = ij.op().run("stats.max", sinusoid)
println("The maximum pixel value:")
println("* before pointwise square root is:\t\t" + imgMax)
println("* after pointwise square root should be:\t" + ij.op().run("math.sqrt", imgMax.get()))
// Create and save a sqrt op operating on DoubleType
sqrtOp = Computers.unary(ij.op(), Ops.Math.Sqrt.class, DoubleType.class, DoubleType.class)
// 'Map' the op on the 'sinusoid' image to calculate the square root of every pixel.
ij.op().run("map", mapOutput, sinusoid, sqrtOp);
println("* after pointwise square root is:\t\t" + ij.op().stats().max(mapOutput))
""
The maximum pixel value: * before pointwise square root is: 19.999118601072674 * after pointwise square root should be: 4.472037410518015 * after pointwise square root is: NaN
In the above example, we additionally calculate some image statistics, namely the maximum pixel value in order to verify that the map
op is actually doing what we'd expect. First we print the maximum of the sinusoid
image and the square root of that value. After the map operation, the new maximum of the map output should be the same as previously calculated one, which is indeed the case :)
Next, we want to look at a slightly more advanced example of chaining. Here, we want to calculate the Difference of Gaussian of an image which is a classical image processing operation. If we look at the method signatures of the dog
op, we see that there is a version with the following inputs:
Therefore, we will create all of these ops and call the dog
op using them:
import net.imglib2.type.numeric.real.DoubleType
import net.imagej.ops.special.computer.Computers
import net.imagej.ops.Ops;
import net.imglib2.RandomAccessibleInterval
import java.util.Collections;
import net.imagej.ops.special.function.Functions;
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory;
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory.Boundary;
// For this example we will use a different image as this will illustrate the result of the `dog` more nicely
input = ij.scifio().datasetIO().open("https://imagej.net/images/clown.png")
// Convert the input to float in order to get a corerct result
input = ij.op().run("convert.float32", input)
Class raiClass = RandomAccessibleInterval.class
// Create two Gaussian filter ops with different sigmas as secondary inputs. Furthermore, we need to specify an out-of-bounds strategy
gaussSigma2Op = Computers.unary(ij.op(), Ops.Filter.Gauss.class, input, input, 2, new OutOfBoundsMirrorFactory(Boundary.SINGLE))
gaussSigma8Op = Computers.unary(ij.op(), Ops.Filter.Gauss.class, input, input, 8, new OutOfBoundsMirrorFactory(Boundary.SINGLE))
// Create the tmp and output creator ops
outputCreatorFunc = Functions.unary(ij.op(), Ops.Create.Img.class, raiClass, input)
tmpCreatorFunc = Functions.unary(ij.op(), Ops.Create.Img.class, raiClass, input)
// Allocate output
dogOutput = ij.op().run("create.img", input)
// Run the `dog` op
ij.op().run("dog", dogOutput, input, gaussSigma2Op, gaussSigma8Op, outputCreatorFunc, tmpCreatorFunc)
ij.notebook().display(dogOutput)