ImgLib2: any way to apply a math operation to three RandomAccessibleInterval and store the result in a fourth, without an explicit loop?

ImgLib2: any way to apply a math operation to three RandomAccessibleInterval and store the result in a fourth, without an explicit loop?

The LoopBuilder seems limited to two inputs only.

If you convert them into IterableInterval’s, then you can use fun.imagej.img/map-img with a lambda that takes 4 IterableInterval’s as arguments. However, under the hood we’re expanding to a Clojure loop, but you might enjoy it :smiley: I am more than happy to provide specific Clojure snippets (and if I can’t provide them off the cuff, then I’ll try to add support as needed).

1 Like

And for the sake of completeness:

(ns albert-loves-clojure
  (:require [fun.imagej.img :as img]
            [fun.imagej.ops]
            [fun.imagej.img.cursor :as cursor]
            [fun.imagej.core :as ij]))

(let [img1 (fun.imagej.ops.create/img (long-array [25 25]))
      img2 (fun.imagej.ops.create/img (long-array [25 25]))
      img3 (fun.imagej.ops.create/img (long-array [25 25]))
      img4 (fun.imagej.ops.create/img (long-array [25 25]))]
  (img/map-img (fn [cur1 cur2 cur3 cur4]
                  (cursor/set-val cur4 (* 255 (java.lang.Math/random))))
                img1 img2 img3 img4)
  (ij/show img4))
3 Likes

Another option, albeit a little cumbersome, would be to merge the three input RAI into a single RAI, either by applying Views.pair twice, or by calling Views.collapse( Views.stack ). You would then have to dismantle the Pair or Composite for your actual function call, of course.

That’s interesting. Care to provide a small example?

@hanlovsky, I’ve pushed a branch that would do what I need, called “pixel-wise”. Care to comment? I’d like to submit it as a pull request, but the button is greyed out.

https://github.com/imglib/imglib2/tree/pixel-wise

See classes LoopMath and LoopMathTest.

https://github.com/imglib/imglib2/blob/pixel-wise/src/main/java/net/imglib2/loops/LoopMath.java
https://github.com/imglib/imglib2/blob/pixel-wise/src/test/java/net/imglib2/loops/LoopMathTest.java

import net.imglib2.img.array.ArrayImg;
import net.imglib2.img.array.ArrayImgs;
import net.imglib2.img.basictypeaccess.array.IntArray;
import net.imglib2.loops.LoopBuilder;
import net.imglib2.type.numeric.integer.IntType;
import net.imglib2.util.Pair;
import net.imglib2.view.Views;
import net.imglib2.view.composite.RealComposite;

import java.util.Arrays;
import java.util.function.BiConsumer;

public class Example {

	static void someFunction( IntType i1, IntType i2, IntType i3, IntType o )
	{
		o.set( i1.get() + i2.get() + i3.get() );
	}

	public static void main( String[] args ) {

		ArrayImg< IntType, IntArray > im1 = ArrayImgs.ints( new int[]{ 0, 1 }, 2 );
		ArrayImg< IntType, IntArray > im2 = ArrayImgs.ints( new int[]{ 1, 2 }, 2 );
		ArrayImg< IntType, IntArray > im3 = ArrayImgs.ints( new int[]{ 2, 3 }, 2 );
		ArrayImg< IntType, IntArray > im4 = ArrayImgs.ints( new int[]{ 0, 0 }, 2 );
		ArrayImg< IntType, IntArray > im5 = ArrayImgs.ints( new int[]{ 0, 0 }, 2 );

		LoopBuilder< LoopBuilder.TriConsumer< IntType, Pair< IntType, IntType >, IntType > > lb1 = LoopBuilder.setImages( im1, Views.interval( Views.pair( im2, im3 ), im2 ), im4 );
		LoopBuilder< BiConsumer< RealComposite< IntType >, IntType > > lb2 = LoopBuilder.setImages( Views.collapseReal( Views.stack( im1, im2, im3 ) ), im5 );

		lb1.forEachPixel( (i1, p, i4) -> someFunction(i1, p.getA(), p.getB(), i4 ) );
		lb2.forEachPixel( (c, i5) -> someFunction(c.get( 0 ), c.get( 1 ), c.get( 2 ), i5 ) );

		System.out.println( Arrays.toString(im4.update( null).getCurrentStorageArray() ) );
		System.out.println( Arrays.toString(im5.update( null).getCurrentStorageArray() ) );

	}

}

Output is

[3, 6]
[3, 6]

Thanks a lot!

To comment that it seems that Views.pair is like python’s zip function, delivering a stream of tuples.

It would be nice if this tuple could be of any size, not just size 2.

I saw you were able to create a PR now, so I am sure that people will comment on that.

The PR aims at addressing the final part: defining the collating function in java itself, rather than in a scripting language.

I am sure it can be done better, this PR is a starting point.
Also, it could be trivially parallelized.

That’s what Views.collapse does, with the restrictions, that all types have to be the same. If you would like to have more freedom in the types that you can pass, you would have to use dynamic casts (bad) or explicitly reproduce the Views.pair behavior for any number of elements (triple, …) individually as Java does not offer variadic generics (like variadic templates in C++).

3 Likes

Hi, just found this here after closing the pull request. Too many forums, chats, repos without automatic cross-references… I think it’s not ready for a pull request (which implies that merging at any time soon is appropriate), and it’d be fine to have it sit around as a topic branch until it feels ready. I also think that building it with multiple loops and loop-builder would be a better way because that may make it more efficient under real world conditions with multiple types and backends mixed (@maarzt, @tpietzsch ?) . I would also appreciate comments from the ops folks because I always believed that this particular scenario is one of the core things that imagej-ops was built for? @dietzc @ctrueden?

1 Like

Hi Stephan,

I’ve pushed further changes to the “pixel-wise” branch to rename classes to avoid confusion with ops, and to implement the Converter so that values from the input image types are passed onto the output in a specifiable way, rather than by setReal(getRealDouble()).

I see this LoopMath class as actually better than the LoopBuilder, which does not address the key looping problem from the scripting perspective: not only the looping itself, but also the execution of the function at every pixel. LoopBuilder can simplify java code. LoopMath aims at being far more general than that: so it is also more limited (6 math operations only).

Question: if you still like this class but don’t want it in the “main” repos, despite the similar LoopBuilder being there, then, in which repository would you accept this LoopMath class?

By the way I’ve tested the performance and, despite the many motions compared to plain looping and math, it’s within 1.5X (as the unit test shows). It is also trivially parallelizable, haven’t implemented it yet.

Also, I’ve commented on the pull request, answering the points you raised. Given that you closed it without awaiting a response, I am not sure that you saw it. The LoopMath won’t run if images differ in dimension number, size, or iterating order, throwing an Exception with that information.

I am convinced that this single class makes a real difference towards being able to use ImgLib2 from a scripting environment, just like the Views, RealViews and Converters were total game changers. See select examples here using the latter:

https://www.ini.uzh.ch/~acardona/fiji-tutorial/#s11

Best,

Albert

1 Like
ij.op().run("math.multiply", image1, image2)

Seems pretty straight forward

Also please note that there is a map op in imagej-ops: https://github.com/imagej/tutorials/blob/op-directory/notebooks/1_-_Using_ImageJ/Ops/map.ipynb

Ops also has a bunch of type matching stuff as well, it would be a shame if you didn’t try it out first.

@albertcardona Without having looked at your code in detail, how is introducing a LoopMath class different from offering specific implementations of Consumer that can be passed into LoopBuilder.forEachPixel?

@kephale : the ij.op() library is powerful, but with this power comes an additional complexity that makes it hard to teach to the casual coder like e.g. biology students (disclosure: I am a biology student myself). It is not as composable either, at least not in a pretty, simple way. LoopMath aims at being simple (only 6 maths operations: the same that are supported under net.imglib2.type.operators, in addition to the NumericType.compareTo for Min and Max), and at being composable.

@hanslovsky : The problem with passing a Consumer created in e.g. jython to the LoopBuilder.forEachPixel is that the jython interpreter is slow, so the execution of the body of the Consumer is slow. Keeping the math entirely within the java side of things is orders of magnitude more performant.

An example of code with LoopMath: compute the brightness of an RGB image:

RandomAccesibleInterval<ARGBType> img = ...
RandomAccesibleInterval<UnsignedByteType> red = Converters.argbChannel(img, 1);
RandomAccesibleInterval<UnsignedByteType> green = Converters.argbChannel(img, 2);
RandomAccesibleInterval<UnsignedByteType> blue = Converters.argbChannel(img, 3);
ArrayImg< FloatType, ? > brightness = new ArrayImgFactory< FloatType >( new FloatType() ).create( img );
LoopMath.compute( brightness, new Div( new Max( red, new Max( green, blue ) ), 3.0 ) );

In Jython:

img = ...
red = Converters.argbChannel(img, 1);
green = Converters.argbChannel(img, 2);
blue = Converters.argbChannel(img, 3);
brightness = ArrayImgFactory( FloatType() ).create( img )
LoopMath.compute( brightness, Div( Max( red, Max( green, blue ) ), 3.0 ) )

In Jython, the above runs within 1.5X the time of a plain unfolded loop in java (see unit tests).

Sorry for the confusion @albertcardona, I was actually referring to providing a set of Consumers that are implemented in Java, thus circumventing the slow interpreter issue.

If composition (potentially with shortcuts) is of importance, that sounds like a separate project, something like imglib2-compute-graph.

As a side note (knowing that this is just example code): Looking at what Converters.argbChannel does, I would advise against using it for read-only operations as the ARGBChannelSamplerConverter creates a new UnsignedByteType for every single pixel.

Thanks @hanslovsky, is imglib2-compute-graph publicly available?

On the Consumers: isn’t in a way what the LoopMath.Add, .Sub, .Mul, .Div, .Max, .Min are? And they are composable with each other and with images and plain scalar numbers.

On the ARGBChannelSamplerConverter: the instantiation of a new UnsignedByteType for every single pixel seems like an issue that could be fixed.

There is no such thing like imglib2-compute-graph, yet. My thought is that, if something like concatenated operations is to be introduced, it should be done in a well-thoughout manner in a separate package for compute graphs.

Consumers: Yes, that’s exactly what they are and I do not see a reason to introduce a new LoopMath framework if we could just have consumers to pass into the LoopBuilder framework. Composition is nice to have but, as stated above, seems like a more involved project if done properly in the context of a compute graph. If you look at CPython and numpy, they do not compose operations either, instead you write something like a = b + c * d and each operation is executed in its own loop.

ARGBChannelSamplerConverter: I don’t think there is a (simple) way around it if you want to have both the read and the write direction (for read-only, it is trivial). A proxy-object pool like in the mastodon-collection could do the trick but adds complexity.

Thanks for the explanations, @hanslovsky

Is there a read-only alternative to ARGBChannelSamplerConverter which avoids creating a new ARGBType for every pixel access?