imglib / imglib2-algorithm

Image processing algorithms for ImgLib2
http://imglib2.net/
Other
22 stars 20 forks source link

ImgMath: new package net.imglib2.algorithm.math #70

Closed acardona closed 6 years ago

acardona commented 6 years ago

This new package contains only one class, ImgMath. Used to be called "LoopMath" in the pixel-wise branch of the main imglib2 repository (see https://github.com/imglib/imglib2/pull/221). Moving it to imglib2-algorithm by suggestion of @axtimwalde.

The core idea is to be able to perform pixel-wise math without an explicit loop, which is expensive in some scripting languages such as Jython.

A possible alternative exist, such as LoopBuilder, but it suffers from having to define the function to execute at every pixel in the scripting language, which can be orders of magnitude slower than in java. In addition, currently LoopBuilder accepts up to 2 input images only, and composability of operations is verbose.

This new package is not a replacement for ops. It is far, far simpler than ops, and limited to 7 operations: Add, Sub, Mul, Div, Min, Max, Neg. These operations cover a surprisingly large amount of relatively simple operations that come up in everyday image processing. The simplicity of ImgMath dispenses with the need to learn the feature-rich ops package for simple image arithmetic.

All math operations are Type-based, using the interfaces implemented by NumericType. Currently, the most abstract Type class that can be handled is RealType.

The input images can have different iteration modes. When so, ImgMath will loop through every pixel using a Localizable approach. Otherwise it will use the Cursor/Iterable approach.

When input images have different number of dimensions, or dimension size doesn't match, an Exception is thrown with an appropriate message.

Top performance is within 1.5X of manually implementing low-level loops with Cursor in java (see unit tests), but the mean run per iteration (over 30 iterations) is the same (~1X), and the max is surprising much lower at 0.25x.

An example, in jython:

from net.imglib2.converter import Converters
from net.imglib2.img.display.imagej import ImageJFunctions as IL
from net.imglib2.img.array import ArrayImgFactory
from net.imglib2.type.numeric.real import FloatType
from net.imglib2.algorithm.math import ImgMath
from net.imglib2.algorithm.math.ImgMath import Max, Div
from ij import IJ

# Load an RGB or ARGB image
imp = IJ.getImage()

# Access its pixel data from an ImgLib2 data structure:
# a RandomAccessibleInterval<ARGBType>
img = IL.wrapRGBA(imp)

# Read out single channels
red = Converters.argbChannel(img, 1)
green = Converters.argbChannel(img, 2)
blue = Converters.argbChannel(img, 3)

# Create an empty image of type FloatType (floating-point values)
# Here, the img is used to read out the interval: the dimensions for the new image
brightness = ArrayImgFactory(FloatType()).create(img)

# Compute the brightness: pick the maximum intensity pixel of every channel
# and then normalize it by dividing by the number of channels
ImgMath(Div(Max([red, green, blue]), 255.0)).into(brightness)

# Show the brightness image
impB = IL.wrap(brightness, imp.getTitle() + " brightness")
impB.show()
acardona commented 6 years ago

The error in travis has nothing to do with the img-math branch. It's related to LocalExtrema.

imagejan commented 6 years ago

@acardona the [ERROR] reporting from maven-javadoc is misleading, as most of these are just warnings. The only actual error is indeed in ImgMath, where an Exception is declared in the javadoc, but not thrown in the code:

/home/travis/build/imglib/imglib2-algorithm/src/main/java/net/imglib2/algorithm/math/ImgMath.java:138: error: exception not thrown: java.lang.Exception
[ERROR]      * @throws Exception When images have different dimensions.

https://github.com/imglib/imglib2-algorithm/blob/c26d083109446b83c53f31684c3eb9e66b329ea5/src/main/java/net/imglib2/algorithm/math/ImgMath.java#L138-L140

acardona commented 6 years ago

Thanks @imagejan you are right. Just fixed the issue.

acardona commented 6 years ago

Here is a minimal example illustrating the use of the "img-math" branch, in jython. The code in java would look nearly identical (except for varags calls, see below).

The computation can be defined with the constructors (Compute, Max, Div) or via the convenient homonimous static methods in the ImgMath class (e.g. compute, max, div) which look more like plain functions and, in java, do not require the cumbersome "new" keyword.

The first compute is simple: compute the brightness by finding the max value of the 3 color channels at each pixel, and dividing by 255 to normalize between 0 and 1.

The second compute, for the image color saturation, showcases a number of functions, including a let for declaration of within-loop variables, and an IF / THEN ELSE statement. Note that THEN and ELSE are identity functions and can be omitted: they are not necessary, and are merely used here for readability. The IF takes 3 arguments, and the THEN (first argument) is executed if the computed value for its first argument is non-zero, with the ELSE (second argument) being executed otherwise. (If you use ELSE first and THEN second, it will work just the same but confuse you to no end, as the THEN would be executed when false, and the ELSE when true).

let declarations can be nested: the most inner one wins. Internally, that's what let does when more than one variable is declared: nest the declarations in multiple class Let instances. So variables (read with var) can be redefined.

Other functions not used here include GT (class GreaterThan), LT (class LessThan), which set the output pixel value to 1 when first arg is larger or smaller than the 2nd, respectively, and to 0 otherwise. EQ (class Equal) sets it to 1 when true, 0 otherwise. And NEQ (class NotEqual) the opposite.

Notice e.g. the max function (class Max) and all other ABinaryFunction subclasses (2-operation functions), including the Let, are capable of varargs, which from jython must be accessed with an array (hence the square brackets in the code below). Otherwise defining max( img1, max( img2, img3 ) ) also works, but it is less readable. (Behind the curtains, that's what the e.g. Max class is doing: instantiating a chain of binary comparisons.)

Compared to using LoopBuilder with a jython-defined function, the performance is excellent: within 1.5-2X of low-level java imglib2 for computations not involving IFstatements, and within 3X for computations with an IF.

from net.imglib2.converter import Converters
from net.imglib2.img.display.imagej import ImageJFunctions as IL
from net.imglib2.img.array import ArrayImgFactory
from net.imglib2.type.numeric.real import FloatType
from net.imglib2.algorithm.math.ImgMath import compute, max, min, div, let, var, IF, THEN, ELSE, EQ, sub
from ij import IJ

# Load an RGB or ARGB image
imp = IJ.openImage("https://imagej.nih.gov/ij/images/leaf.jpg") 
#imp = IJ.getImage()

# Access its pixel data from an ImgLib2 data structure:
# a RandomAccessibleInterval<ARGBType>
img = IL.wrapRGBA(imp)

# Read out single channels
red = Converters.argbChannel(img, 1)
green = Converters.argbChannel(img, 2)
blue = Converters.argbChannel(img, 3)

# Create an empty image of type FloatType (floating-point values)
# Here, the img is used to read out the interval: the dimensions for the new image
brightness = ArrayImgFactory(FloatType()).create(img)

# Compute the brightness: pick the maximum intensity pixel of every channel
# and then normalize it by dividing by the number of channels
compute(div(max([red, green, blue]), 255.0)).into(brightness)

# Show the brightness image
impB = IL.wrap(brightness, imp.getTitle() + " brightness")
impB.show()

# Compute now the image color saturation
saturation = ArrayImgFactory(FloatType()).create(img)
compute( let("red", red,
                     "green", green,
                     "blue", blue,
                     "max", max([var("red"), var("green"), var("blue")]),
                     "min", min([var("red"), var("green"), var("blue")]),
                     IF ( EQ( 0, var("max") ),
                       THEN( 0 ),
                       ELSE( div( sub(var("max"), var("min")),
                                         var("max") ) ) ) ) ).into(saturation)

impB = IL.wrap(saturation, imp.getTitle() + " saturation")
impB.show()
acardona commented 6 years ago

Eclipse and javac do not complain about this "error":

[ERROR] COMPILATION ERROR : [INFO] ------------------------------------------------------------- [ERROR] /home/travis/build/imglib/imglib2-algorithm/src/main/java/net/imglib2/algorithm/math/Compute.java:[84,112] incompatible types: net.imglib2.Cursor cannot be converted to net.imglib2.Cursor [INFO] 1 error

But Travis does. All changes are committed.

The "error" is a cast that travis doesn't like.

ctrueden commented 6 years ago

That error also happens on my macOS system from the CLI via Maven. It is a real error that needs to be fixed or worked around.

acardona commented 6 years ago

The error in travis was fixed by e0c39a1, which rewrites the type constraints in the failing function. Essentially, nothing changed at all in the code; this is a work around a compiler issue in travis.

ctrueden commented 6 years ago

It is important to never use RealType<?> or ? extends RealType<?>. Everywhere in the codebase that does so has substantial typing issues, because it breaks the recursive typing constraint. Even use of raw types works only limitedly, depending on the situation, when the generic bounds are recursive.

Quickly glancing at the code, my impression is that to fix this issue properly, generic bounds need to be added to Compute and IFunction and probably many other places that all use RealType<?>.

acardona commented 6 years ago

The goal of the library is to avoid explicit loops or in-loop function definitions, and to be easy to read and write for humans while remaining high performance when executed from scripting languages or java. All these goals have been met.

Specifying the types is not even possible without a Converter per input image, which is unnecessary as the only guarantee necessary is the RealType interface. And a Converter per input type would be a significant performance drag.

The consumers of this library need not worry about type safety at all beyond the RealType constraint.

The library is feature-complete as is and ready to merge.

On Aug 1, 2018, at 5:51 PM, Curtis Rueden notifications@github.com wrote:

It is important to never use RealType<?> or ? extends RealType<?>. Everywhere in the codebase that does so has substantial typing issues, because it breaks the recursive typing constraint. Even use of raw types works only limitedly, depending on the situation, when the generic bounds are recursive.

Quickly glancing at the code, my impression is that to fix this issue properly, generic bounds need to be added to Compute and IFunction and probably many other places that all use RealType<?>.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

ctrueden commented 6 years ago

The consumers of this library need not worry about type safety at all beyond the RealType constraint.

Then I would suggest you use the raw type RealType instead of RealType<?> or ? extends RealType<?>.

@tpietzsch, what do you think?

acardona commented 6 years ago

OK to merge? @axtimwalde @tpietzsch @StephanPreibisch ?

Can do the merge myself if there aren't any objections.

axtimwalde commented 6 years ago

I am fine with that.