Display an image with non-linear calibrated axes

calibrated-axes-transform.groovy

#@ DatasetService ds
#@ net.imagej.units.UnitService units
#@ UIService ui

// Parameters to play with!

// log-linear equation: cal(raw) = a + b * ln(c + d*raw)
xa = 0; xb = 10; xc = 1; xd = 1   // cal(xRaw) = 10*ln(1+xRaw)
ya = 0; yb = 5000; yc = 1; yd = 3000 // cal(yRaw) = 5000*ln(1+3000*yRaw)
// units
xUnit = "micron"
yUnit = "nm"
// image dimensions
w = 10 // [0, 9] -> [0, ~22] microns
h = 15 // [0, 14] -> [0, ~53227] nm -> [0, ~53] microns

// Construct a plain uncalibrated image
import net.imglib2.img.array.ArrayImgs
data = new byte[w*h];
for (int i=0; i<data.length; i++) data[i] = i;
img = ArrayImgs.unsignedBytes(data, w, h)

import net.imagej.DefaultDataset;
import net.imagej.ImgPlus;
import net.imagej.axis.Axes;
import net.imagej.axis.CalibratedAxis;
import net.imagej.axis.LogLinearAxis;

String name = "Raw log-scale data";
CalibratedAxis xAxis = new LogLinearAxis(Axes.X, xUnit, xa, xb, xc, xd);
CalibratedAxis yAxis = new LogLinearAxis(Axes.Y, yUnit, ya, yb, yc, yd)
ImgPlus imgPlus = new ImgPlus(img, name, xAxis, yAxis)

ui.show(imgPlus);

// RandomAccessibleInterval -> RealRandomAccessible
import net.imglib2.view.Views;
import net.imglib2.RealRandomAccessible;
import net.imglib2.interpolation.randomaccess.NearestNeighborInterpolatorFactory;
interpolator = new NearestNeighborInterpolatorFactory();
RealRandomAccessible realRaw = Views.interpolate(imgPlus, interpolator);

// Implement an invertible RealTransform that respects the axis calibrations and units.

import net.imagej.units.UnitService;
import net.imglib2.RealLocalizable;
import net.imglib2.RealPositionable;
import net.imglib2.realtransform.InverseRealTransform;
import net.imglib2.realtransform.InvertibleRealTransform;

public class CalibratedAxesTransform implements InvertibleRealTransform {

	private UnitService units;

	/** Common destination unit for transformed coordinates. */
	private String unit;

	private CalibratedAxis[] axes;

	public CalibratedAxesTransform(UnitService units, CalibratedAxis... axes) {
		this(units, axes[0].unit(), axes);
	}

	public CalibratedAxesTransform(UnitService units, String unit, CalibratedAxis... axes) {
		this.units = units;
		this.unit = unit;
		this.axes = axes;
	}

	@Override
	public int numSourceDimensions() { return axes.length; }

	@Override
	public int numTargetDimensions() { return axes.length; }

	/** For use with {@link #apply}. */
	private double calibrated(final CalibratedAxis axis, final double raw) {
		// raw -> calibrated
		double calibrated = axis.calibratedValue(raw);

		// source unit -> target unit
		if (unit != null) calibrated = units.value(calibrated, axis.unit(), unit);

		return calibrated;
	}

	/** For use with {@link #applyInverse}. */
	private double raw(final CalibratedAxis axis, double calibrated) {
		// target unit -> source unit
		if (unit != null) calibrated = units.value(calibrated, unit, axis.unit());

		// calibrated -> raw
		return axis.rawValue(calibrated);
	}

	@Override
	public void apply(final double[] source, final double[] target) {
		for (int d=0; d<axes.length; d++) {
			target[d] = calibrated(axes[d], source[d]);
		}
	}

	@Override
	public void apply(final RealLocalizable source, final RealPositionable target) {
		for (int d=0; d<axes.length; d++) {
			double raw = source.getDoublePosition(d)
			double calibrated = calibrated(axes[d], raw);
			target.setPosition(calibrated, d);
		}
	}

	@Override
	public void applyInverse(final double[] source, final double[] target) {
		for (int d=0; d<axes.length; d++) {
			source[d] = raw(axes[d], target[d]);
		}
	}

	@Override
	public void applyInverse(final RealPositionable source, final RealLocalizable target) {
		for (int d=0; d<axes.length; d++) {
			double calibrated = target.getDoublePosition(d);
			double raw = raw(axes[d], calibrated);
			source.setPosition(raw, d);
		}
	}

	@Override
	public InvertibleRealTransform inverse() {
		return new InverseRealTransform(this);
	}

	@Override
	public CalibratedAxesTransform copy() {
		return new CalibratedAxesTransform(units, unit, axes);
	}
}

// Whew! Now we can use our transform:

// Construct the transform.
axes = new CalibratedAxis[imgPlus.numDimensions()];
imgPlus.axes(axes);
axesTransform = new CalibratedAxesTransform(units, axes);

// Create a transformed version of the raw image.
import net.imglib2.realtransform.RealViews;
RealRandomAccessible realCalibrated = RealViews.transformReal(realRaw, axesTransform);

// Rasterize our calibrated RRA, putting it back on an integer grid.
import net.imglib2.RandomAccessible;
RandomAccessible integerCalibrated = Views.raster(realCalibrated);

// Bound our calibrated RA, with the same corners as the actual data.
min = new long[axes.length];
max = new long[axes.length];
double[] sourceMin = new double[axes.length];
double[] sourceMax = new double[axes.length];
for (int d=0; d<axes.length; d++) {
	sourceMin[d] = imgPlus.min(d);
	sourceMax[d] = imgPlus.max(d);
}
double[] targetMin = new double[axes.length];
double[] targetMax = new double[axes.length];
axesTransform.apply(sourceMin, targetMin);
axesTransform.apply(sourceMax, targetMax);

long[] min = new long[axes.length];
long[] max = new long[axes.length];
for (int d=0; d<axes.length; d++) {
	min[d] = Math.ceil(targetMin[d]);
	max[d] = Math.floor(targetMax[d]);
}

calibratedRAI = Views.interval(integerCalibrated, min, max);

ui.show("Calibrated RAI", calibratedRAI);