Using ImgLib2 to shear an image

I wanted to reply to this old forum post about StackReg with an ImgLib2-based solution for image shearing.

Here is my first cut:

// @DatasetService datasetService
// @Img image
// @OUTPUT Dataset result

import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory
import net.imglib2.realtransform.AffineTransform2D
import net.imglib2.realtransform.RealViews
import net.imglib2.view.Views

// declare how we want the image to be interpolated
field = Views.interpolate(image, new NLinearInterpolatorFactory())

// shear the image
affine = new AffineTransform2D()
affine.set([[0, 1, 2], [2, 3, 4], [0, 0, 1]] as double[][])
sheared = RealViews.affine(field, affine)

// apply the original bounding box to the sheared image
bounded = Views.interval(sheared, image)

// pad with 0 (for indices with no defined value, due to the shearing)
extended = Views.extendZero(bounded)

// restore the original bounding box again
interval = Views.interval(extended, image)

// convert the result to an ImageJ dataset
result = datasetService.create(interval)

However, when run against the blobs image, this script fails with:

java.lang.ArrayIndexOutOfBoundsException: -767
    at net.imglib2.img.basictypeaccess.array.ByteArray.getValue(ByteArray.java:63)
    at net.imglib2.type.numeric.integer.GenericByteType.getValue(GenericByteType.java:97)
    at net.imglib2.type.numeric.integer.GenericByteType.set(GenericByteType.java:170)
    at net.imglib2.type.numeric.integer.GenericByteType.set(GenericByteType.java:50)
    at net.imglib2.interpolation.randomaccess.NLinearInterpolator2D.get(NLinearInterpolator2D.java:90)
    at net.imglib2.interpolation.randomaccess.NLinearInterpolator2D.get(NLinearInterpolator2D.java:48)
    at net.imglib2.realtransform.AffineRandomAccessible$AffineRandomAccess.get(AffineRandomAccessible.java:180)
    at net.imglib2.util.Util.getTypeFromInterval(Util.java:759)
    at net.imglib2.view.Views.extendZero(Views.java:198)

Of course, I realize my transformation matrix is nonsense. This was just my initial attempt to get any affine transformation working from within ImageJ.

Anybody have any thoughts on what I’m doing wrong? @axtimwalde @tpietzsch @bogovicj ?

I guess the Views.extendZero doesn’t actually “fill in” missing pixels inside the interval bounds. But I don’t know how to do that? And I am also wondering if there is a general way to get the “best bounding box” for the new sheared data, such that it is does not crop out any of the old pixels?

1 Like

It seems to me that the order of your calls to (Real)Views methods is wrong. The source image should be extended and interpolated before applying the transform (RealViews.affine), so the RandomAccess of the transformed image can sample from out of bounds pixels. I wrote a small working example in Java that you can modify according to your needs.

import java.util.Random;

import net.imglib2.RealRandomAccessible;
import net.imglib2.img.array.ArrayImg;
import net.imglib2.img.array.ArrayImgs;
import net.imglib2.img.basictypeaccess.array.DoubleArray;
import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory;
import net.imglib2.realtransform.AffineGet;
import net.imglib2.realtransform.AffineRandomAccessible;
import net.imglib2.realtransform.AffineTransform2D;
import net.imglib2.realtransform.RealViews;
import net.imglib2.type.numeric.real.DoubleType;
import net.imglib2.view.IntervalView;
import net.imglib2.view.Views;


public class WorkingExampleTransform
{

	public static void main( final String[] args )
	{
		final Random rng = new Random( 100 );
		final ArrayImg< DoubleType, DoubleArray > image = ArrayImgs.doubles( 3, 10 );
		for ( final DoubleType i : image )
		{
			i.set( rng.nextDouble() );
		}

		// declare how we want the image to be interpolated
		final RealRandomAccessible< DoubleType > field = Views.interpolate( Views.extendZero( image ), new NLinearInterpolatorFactory<>());

		// shear the image
		final AffineTransform2D affine = new AffineTransform2D() ;
		affine.set(new double[][] {{0., 1., 2.}, {2., 3., 4.}, {0., 0., 1.} });
		final AffineRandomAccessible< DoubleType, AffineGet > sheared = RealViews.affine(field, affine);

		// apply the original bounding box to the sheared image
		final IntervalView< DoubleType > bounded = Views.interval(sheared, image);

		for ( final DoubleType b : bounded )
		{
			System.out.println( b );
		}

	}

}
1 Like

Doh! Thanks for catching that.

Based on your changes, here is my updated working Groovy script:

// @DatasetService datasetService
// @Img image
// @double(label="Shearing factor", value=1) factor
// @OUTPUT Dataset result

import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory
import net.imglib2.realtransform.AffineTransform2D
import net.imglib2.realtransform.RealViews
import net.imglib2.view.Views

// extend the image with zeroes
extended = Views.extendZero(image)

// decide how we want the image to be interpolated
field = Views.interpolate(extended, new NLinearInterpolatorFactory())

// shear the image
affine = new AffineTransform2D()
affine.set([[1, factor, 0], [0, 1, 0], [0, 0, 1]] as double[][])
sheared = RealViews.affine(field, affine)

// apply the original bounding box to the sheared image
bounded = Views.interval(sheared, image)

// convert the result to an ImageJ dataset
result = datasetService.create(bounded)

I would still like to know if there is a super-convenient way to rebound the result to hug the sheared image without cropping it. The script could run the corner coordinates through the affine transform, but that seems like a lot of work—is there a utility method somewhere?

The only thing that I could think of is the BoundingBox class but considers only integer coordinates (not real coordinates), and neither Transform nor RealTransform has a method to transform a bounding box. I agree that a utility method would be useful. However, that will be non-trivial for arbitrary transformations.

Of course. But affines are super common. I will look into creating such a utility method. Well, I’ll probably have one my students do it. But it sounds like it would be useful. :smile:

I filed an issue for it.

1 Like

I had a closer look at the implementation of ShearTransform and found out that it implements BoundingBoxTransform. This interface as it is right now makes sense only for integer Transform. However, it could be extended for RealTransform. The AffineTransform any other transform for which the BoundingBox can be transformed trivially, could then implement this interface.
@axtimwalde and I thought about implementing a real value shear transform with arbitrary shearing factors but decided not to because for our use case we only needed a 45 degree integer shear, which is implemented in ShearTransform. Such a shear can be implemented efficiently with just a single addition.

I will add a comment to the issue, pointing to the BoundingBox class and BoundingBoxTransform interface.

2 Likes

Sorry to invoke this topic after a good long while,
If I apply a 2D transform to a hyper stack (like the example you gave above), the shear would be a 2 dim interval. The conclusion is that the transformation is not applied to every plane then?
In that case, should I manually iterate over every plane of the image ?

The answer would be to define the affine transform n-dimensional!

@Masoudas are you using bigdataviewer ? You can directly define an AffineTransform3D for each source. The affinetransform can unshear the source.

But we need a bit more context in order to give you a useful answer…

2 Likes

Thanks for the response.

I actually need to apply a 2D affine transform to an n-dim Img interface (of ImgLib2). That would imply transforming each plane of the image.

I think extending the 2D transform to an n-dim affine works, right?

Hi @Masoudas,

Again, I’m not completely sure what you want to achieve. If you have some sample code that would probably be easier to guide you. I do not understand the context where you would need to extend the 2D transform.

You have different options in order to perform shearing, but a really important factor is how you want to display the result:

  • Either it is shown in bigdataviewer : then it’s great, because the image is “rasterized” live depending on how the user selects its slicing plane.

Gist: https://gist.github.com/NicoKiaru/f88aaf6bfc201d468b57c771424f1f06

  • or it is exported as a Dataset and shown using IJ1 display, as in the scripts discussed above : here the issue is that you need to rasterize your data yourself, and select the appropriate interval where you expect you data to be. But there’s no need to ‘extend’ the 2D affine transform. The AffineTransfrom3D already exists. I’ve modified the script to make it work in 3D. You can open the sample image “MRI Stack (528K)” to see how it sheared the image in 3D (so all planes).
// @DatasetService datasetService
// @Img image
// @double(label="Shearing factor", value=1) factor
// @OUTPUT Dataset result

import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory
import net.imglib2.realtransform.AffineTransform3D
import net.imglib2.realtransform.RealViews
import net.imglib2.view.Views

// extend the image with zeroes
extended = Views.extendZero(image)

// decide how we want the image to be interpolated
field = Views.interpolate(extended, new NLinearInterpolatorFactory())

// shear the image
affine = new AffineTransform3D()
affine.set([[1, factor, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]] as double[][])
sheared = RealViews.affine(field, affine)

// apply the original bounding box to the sheared image
bounded = Views.interval(sheared, image)

// convert the result to an ImageJ dataset
result = datasetService.create(bounded)


Hope that helps

2 Likes

Thanks for the helpful and very detailed response.
If you check the following code, you see that the idea is to apply a transform to each plane individually, and not just for visualization.

Let me know what you think, or if there’s a better approach (note that I just iterate over the cursor to show that the code works. We better replace the image with an actual 4D image to properly see the result. And I’m looking for applying a generic Affine, not just shear).

 // Create a 4Dim image
    Img<UnsignedShortType> image = new ArrayImgFactory<UnsignedShortType>(new UnsignedShortType()).create(4, 4, 4,
            4);

    final RealRandomAccessible<UnsignedShortType> field = Views.interpolate(Views.extendZero(image),
            new NLinearInterpolatorFactory<>());

    // A 4Dim transformation, that transforms each plane.
    final AffineTransform affine = new AffineTransform(4);
    affine.set(0, 0, 0);
    affine.set(2, 0, 1);
    affine.set(0, 0, 4);

    affine.set(0, 1, 0);
    affine.set(0, 1, 1);
    affine.set(0, 1, 4);

    final AffineRandomAccessible<UnsignedShortType, AffineGet> sheared = RealViews.affine(field, affine);

    // apply the original bounding box to the transformed image
    final IntervalView<UnsignedShortType> bounded = Views.interval(sheared, image);

    // Just printing the position cursor to show that transformation is done properly.
    long[] position = new long[4];
    for (Cursor<UnsignedShortType> cursor = bounded.cursor(); cursor.hasNext();) {
        cursor.next();
        cursor.localize(position);
        System.out.println(position[0] + " " + position[1] + " " + position[2]);
    }
1 Like