Searching for a 3D binary interpolation plugin

plugin

#1

I am searching for a plugin to interpolate between two 3D masks (two 3D binary images). That is, code similar to this by Johannes Schindelin in the VIB library:

… but instead of interpolating in 2D, interpolating in 3D.

The goal is to take 2 arbitrary 3D meshes (where the vertices do not correspond) and produce a sequence of weighted intermediate meshes.

Does it exist? Thanks.


#2

To further add that the BinaryInterpolator class does not seem to be in the same vib package, and I don’t know how to find it. A version of it, adapted to work with ImgLib1 (ancient history!) exists in the TrakEM2 repository as class BinaryInterpolation2D.


#3

@hanslovsky
You’ve done something like this right?


#4

I am not sure that what I did is related:

  • take two 2D slices of label data as input
  • interpolate between those 2D slices (orthogonal to image plane):
    1. signed distance transform for all contained labels on both slices
    2. interpolate each of these distance transforms
    3. take arg min of distance transform as interpolated label

If that is what you need, you will find the code here, and a usage example here and here (the first example is probably going to be more useful for you).


#5

Thanks, it is related, but what I wanted is to interpolate in 3D: given two 3D objects (each represented as a binary mask in 3D), create an interpolated version.

Turns out, it wasn’t so hard. I followed the documentation that I wrote in TrakEM2’s ini.trakem2.imaging.BinaryInterpolation2D class, and wrote one in jython in a few minutes.

If anybody would like to write a java-based implementation, please do. I’d use it. It could go into imglib2-algorithm.

# Albert Cardona 2018-10-22
# A method to generate interpolated masks between two 3D masks
# Should work with any number of dimensions.
# Based on the documentation found in class ini.trakem2.imaging.BinaryInterpolation2D
#
# Note that a java-based implementation would be significantly faster.

from net.imglib2.img.array import ArrayImgs
from org.scijava.vecmath import Point3f
from jarray import zeros, array
from net.imglib2.img.display.imagej import ImageJFunctions as IL
from net.imglib2.view import Views
from net.imglib2 import KDTree, RealPoint
from itertools import imap
from functools import partial
import operator
from net.imglib2.type.numeric.integer import UnsignedByteType
from net.imglib2.neighborsearch import NearestNeighborSearchOnKDTree
from net.imglib2.util import Intervals


# First 3D mask: a sphere
img1 = ArrayImgs.unsignedBytes([100, 100, 100])
p = zeros(3, 'l')
cursor = img1.cursor()
middle = Point3f(49.5,49.5, 49.5)
distance_sq = float(30 * 30)

while cursor.hasNext():
  cursor.fwd()
  cursor.localize(p)
  if middle.distanceSquared(Point3f(p[0], p[1], p[2])) < distance_sq:
    cursor.get().setOne()
  else:
    cursor.get().setZero()

imp1 = IL.wrap(img1, "sphere")
imp1.setDisplayRange(0, 1)
imp1.show()


# Second 3D mask: a cube
img2 = ArrayImgs.unsignedBytes([100, 100, 100])
for t in Views.interval(img2, [20, 20, 20], [80, 80, 80]):
  t.setOne()

imp2 = IL.wrap(img2, "cube")
imp2.setDisplayRange(0, 1)
imp2.show()

# Find edges
def findEdgePixels(img):
  edge_pix = []
  zero = img.firstElement().createVariable()
  zero.setZero()
  imgE = Views.extendValue(img, zero)
  pos = zeros(img.numDimensions(), 'l')
  inc = partial(operator.add, 1)
  dec = partial(operator.add, -1)
  cursor = img.cursor()
  while cursor.hasNext():
    t = cursor.next()
    if 0 == t.getIntegerLong():
      continue
    # Sum neighbors of non-zero pixel: if any is zero, sum is less than 27
    # and we have found an edge pixel
    cursor.localize(pos)
    minimum = map(dec, pos)
    maximum = map(inc, pos) 
    box = Views.interval(imgE, minimum, maximum)
    if sum(imap(UnsignedByteType.getIntegerLong, box)) < 27:
      edge_pix.append(RealPoint(array(list(pos), 'f')))
  return edge_pix

# Generate interpolated image
def makeInterpolatedImage(img1, search1, img2, search2, weight):
  """ weight: float between 0 and 1 """
  img3 = ArrayImgs.unsignedBytes(Intervals.dimensionsAsLongArray(img1))
  c1 = img1.cursor()
  c2 = img2.cursor()
  c3 = img3.cursor()
  while c3.hasNext():
    t1 = c1.next()
    t2 = c2.next()
    t3 = c3.next()
    sign1 = -1 if 0 == t1.get() else 1
    sign2 = -1 if 0 == t2.get() else 1
    search1.search(c1)
    search2.search(c2)
    value1 = sign1 * search1.getDistance() * (1 - weight)
    value2 = sign2 * search2.getDistance() * weight
    if value1 + value2 > 0:
      t3.setOne()
  return img3

edge_pix1 = findEdgePixels(img1)
kdtree1 = KDTree(edge_pix1, edge_pix1)
search1 = NearestNeighborSearchOnKDTree(kdtree1)
edge_pix2 = findEdgePixels(img2)
kdtree2 = KDTree(edge_pix2, edge_pix2)
search2 = NearestNeighborSearchOnKDTree(kdtree2)

steps = []
for weight in [0.2, 0.4, 0.6, 0.8]:
  step = makeInterpolatedImage(img1, search1, img2, search2, weight)
  steps.append(step)

img3 = Views.stack([img1] + steps + [img2])

imp3 = IL.wrap(img3, "interpolated steps")
imp3.setDisplayRange(0, 1)
imp3.show()

(Update: edited to illustrate how to generate a sequence of interpolated steps)


#6

For completeness, here is the documentation from ini.trakem2.imaging.BinaryInterpolation2D:

/** Given two binary images of the same dimensions,
 * generate an interpolated image that sits somewhere
 * in between, as specified by the weight.
 * 
 * For each binary image, the edges are found
 * and then each pixel is assigned a distance to the nearest edge.
 * Inside, distance values are positive; outside, negative.
 * Then both processed images are compared, and wherever
 * the weighted sum is larger than zero, the result image
 * gets a pixel set to true (or white, meaning inside).
 * 
 * A weight of zero means that the first image is not present at all
 * in the interpolated image;
 * a weight of one means that the first image is present exclusively.
 * 
 * The code was originally created by Johannes Schindelin
 * in the VIB's vib.BinaryInterpolator class, for ij.ImagePlus.
 */

#7

Ah, then I misunderstood initially, I thought you were trying to interpolate meshes. 3D should be possible with my code, too, as it is written for arbitrary dimensions but I only tested in 2D. My code is for an arbitrary number of objects and a single object is just a special case of that. I see that you have your own implementation already, though.

label-utilities is available through the imagej maven repo


#8

Thanks for the pointers @hanslovsky. The script above can also handle multiple objects simultaneously. The IDT approach is quite powerful.


#9

Nice that you have this done for n-dimensions @albertcardona - it would be a handy thing to have as an Op in Java. Once the binary masks are calculated one can easily make meshes from them.

Some time ago we used Johannes’ signed distance transform code in BoneJ as Interpolate ROIs , which was more-or-less a bundled version of VIB-lib/vib/BinaryInterpolator.java accessed by a userland plugin that interpolated along the z-axis between 2D ROIs drawn on different stack slices stored in the ROI Manager. Wayne upstreamed it to ImageJ1 and so it has been for several years a component of the ROI Manager.


#10

Thanks @mdoube, that’s great to know that Wayne Rasband upstreamed Johannes Schindelin’s code into ImageJ1. Now I understand why the vib library doesn’t have it anymore.

Upon looking at @hanslovsky code examples, I think his library does the same thing, and multithreaded and with more options. Perhaps that’s the better target for an easier, more accessible library entry point that would act as a wrapper over his code.

If you want the above jython script in java, I am sure it is easy to translate. For performance reasons I rewrote the hot bits in java, below, since I wanted to run it in images that are far larger than the examples shown in the script. Feel free to turn this code into pure java and upload it somewhere where it can be used within Fiji:

# Albert Cardona 2018-10-22
# A method to generate interpolated masks between two 3D masks
# Should work with any number of dimensions.
# Based on the documentation found in class ini.trakem2.imaging.BinaryInterpolation2D
#
# Given two binary images of the same dimensions,
# generate an interpolated image that sits somewhere
# in between, as specified by the weight.
# 
# For each binary image, the edges are found
# and then each pixel is assigned a distance to the nearest edge.
# Inside, distance values are positive; outside, negative.
# Then both processed images are compared, and wherever
# the weighted sum is larger than zero, the result image
# gets a pixel set to true (or white, meaning inside).
# 
# A weight of zero means that the first image is not present at all
# in the interpolated image;
# a weight of one means that the first image is present exclusively.
# 
# The code was originally created by Johannes Schindelin
# in the VIB's vib.BinaryInterpolator class, for ij.ImagePlus,
# but the code below uses a completely different implementation
# strategy based on a KDTree for finding the edges of the masks.
#
#
# Note that a java-based implementation would be significantly faster.

from net.imglib2.img.array import ArrayImgs
from org.scijava.vecmath import Point3f
from jarray import zeros, array
from net.imglib2.img.display.imagej import ImageJFunctions as IL
from net.imglib2.view import Views
from net.imglib2 import KDTree, RealPoint, Cursor, RandomAccessible
from itertools import imap
from functools import partial
import operator
from net.imglib2.type.numeric.integer import UnsignedByteType
from net.imglib2.neighborsearch import NearestNeighborSearchOnKDTree
from net.imglib2.util import Intervals
from fiji.scripting import Weaver
from java.util import ArrayList
from net.imglib2.img import Img


# First 3D mask: a sphere
img1 = ArrayImgs.unsignedBytes([100, 100, 100])
p = zeros(3, 'l')
cursor = img1.cursor()
middle = Point3f(49.5,49.5, 49.5)
distance_sq = float(30 * 30)

w = Weaver.method("""
public void fillSphere(final Cursor cursor, final Point3f middle, final float distance_sq) {
  final long[] p = new long[cursor.numDimensions()];
  while (cursor.hasNext()) {
    cursor.fwd();
    cursor.localize(p);
    final UnsignedByteType t = (UnsignedByteType) cursor.get();
    if (middle.distanceSquared(new Point3f( (float)p[0], (float)p[1], (float)p[2] ) ) < distance_sq) {
      t.setOne();
    } else {
      t.setZero();
    }
  }
}

public void fillOnes(final Cursor cursor) {
  while (cursor.hasNext()) {
    ((UnsignedByteType)cursor.next()).setOne();
  }
}
""", [Cursor, Point3f, UnsignedByteType])

w.fillSphere(img1.cursor(), middle, distance_sq)

imp1 = IL.wrap(img1, "sphere")
imp1.setDisplayRange(0, 1)
imp1.show()


# Second 3D mask: three small cubes
img2 = ArrayImgs.unsignedBytes([100, 100, 100])

w.fillOnes(Views.interval(img2, [10, 10, 10], [29, 29, 29]).cursor())
w.fillOnes(Views.interval(img2, [70, 10, 70], [89, 29, 89]).cursor())
w.fillOnes(Views.interval(img2, [40, 70, 40], [59, 89, 59]).cursor())

imp2 = IL.wrap(img2, "cube")
imp2.setDisplayRange(0, 1)
imp2.show()

w2 = Weaver.method("""
public final ArrayList findEdgePixels(final Img img) {
  final ArrayList edge_pix = new ArrayList();
  final RandomAccessible imgE = Views.extendValue(img, ((UnsignedByteType)img.firstElement()).createVariable());
  final long[] pos = new long[img.numDimensions()];
  final long[] minimum = new long[pos.length],
               maximum = new long[pos.length];
  Cursor cursor = img.cursor();
  while (cursor.hasNext()) {
    final UnsignedByteType t = (UnsignedByteType) cursor.next();
    if (0 == t.getIntegerLong()) continue;
    cursor.localize(pos);
    for (int i=0; i<pos.length; ++i) minimum[i] = pos[i] - 1;
    for (int i=0; i<pos.length; ++i) maximum[i] = pos[i] + 1;
    final Cursor ce =  Views.interval(imgE, minimum, maximum).cursor();
    while (ce.hasNext()) {
      final UnsignedByteType tv = (UnsignedByteType) ce.next();
      if ( 0 == tv.getIntegerLong() ) {
        final float[] fp = new float[pos.length];
        for (int i=0; i<pos.length; ++i) fp[i] = (float) pos[i];
        edge_pix.add(new RealPoint(fp));
        break;
      }
    }
  }
  return edge_pix;
}

public final Img makeInterpolatedImage(final Img img1, final NearestNeighborSearchOnKDTree search1,
                                       final Img img2, final NearestNeighborSearchOnKDTree search2,
                                       final float weight) {
  final Img img3 = ArrayImgs.unsignedBytes(Intervals.dimensionsAsLongArray(img1));
  final Cursor c1 = img1.cursor(),
               c2 = img2.cursor(),
               c3 = img3.cursor();
  while (c3.hasNext()) {
    final UnsignedByteType t1 = (UnsignedByteType) c1.next(),
                           t2 = (UnsignedByteType) c2.next(),
                           t3 = (UnsignedByteType) c3.next();
    final double sign1 = (0 == t1.get() ? -1 : 1),
                 sign2 = (0 == t2.get() ? -1 : 1);
    search1.search(c1);
    search2.search(c2);
    final double value1 = sign1 * search1.getDistance() * (1 - weight),
                 value2 = sign2 * search2.getDistance() * weight;
    if (value1 + value2 > 0) {
      t3.setOne();
    }
  }
  return img3;
}

""", [Img, Cursor, ArrayList, UnsignedByteType,
      Views, RandomAccessible, RealPoint,
      NearestNeighborSearchOnKDTree, ArrayImgs,
      Intervals])

edge_pix1 = w2.findEdgePixels(img1)
kdtree1 = KDTree(edge_pix1, edge_pix1)
search1 = NearestNeighborSearchOnKDTree(kdtree1)
edge_pix2 = w2.findEdgePixels(img2)
kdtree2 = KDTree(edge_pix2, edge_pix2)
search2 = NearestNeighborSearchOnKDTree(kdtree2)

steps = []
for weight in [x / 10.0 for x in xrange(2, 10, 2)]:
  steps.append(w2.makeInterpolatedImage(img1, search1, img2, search2, weight))

imp3 = IL.wrap(Views.stack([img1] + steps + [img2]), "interpolations")
imp3.setDisplayRange(0, 1)
imp3.show()

The reason the java code above reads unindiomatic lays in the limitations of the java compiler in the Tools.jar, which doesn’t know about generics.


#11

Thanks for the praise, I have to lower the expectations a little bit, though: The method itself is not parallelized, only one of the usage examples parallelizes over the pairs of sections between which to interpolate. This should work in 3D, as well, but I haven’t tried so far.

As an example, this is the output of the first usage example:

  • top: first section
  • bottom: second section
  • middle: interpolation with varying weights

interpolated-stack-with-first-and-last