diff --git a/render-app/pom.xml b/render-app/pom.xml index 73bb340b3..fb17b1d4f 100644 --- a/render-app/pom.xml +++ b/render-app/pom.xml @@ -164,6 +164,28 @@ n5-hdf5 + + + org.janelia.saalfeldlab + n5-universe + + + + + org.janelia + n5-ng-precomputed + 0.1.0 + + + + + org.janelia + n5-ng-precomputed + 0.1.0 + tests + test + + diff --git a/render-app/src/main/java/org/janelia/alignment/loader/N5SliceLoader.java b/render-app/src/main/java/org/janelia/alignment/loader/N5SliceLoader.java index d2801f456..cd176707d 100644 --- a/render-app/src/main/java/org/janelia/alignment/loader/N5SliceLoader.java +++ b/render-app/src/main/java/org/janelia/alignment/loader/N5SliceLoader.java @@ -9,9 +9,11 @@ import java.net.URI; import java.net.URLDecoder; import java.nio.charset.Charset; +import java.util.Optional; import net.imglib2.loops.LoopBuilder; import net.imglib2.type.numeric.integer.UnsignedShortType; +import org.janelia.alignment.util.QueryKeyValueParameters; import org.janelia.saalfeldlab.n5.DataType; import org.janelia.saalfeldlab.n5.DatasetAttributes; import org.janelia.saalfeldlab.n5.N5FSReader; @@ -67,44 +69,19 @@ public ImageProcessor load(final String urlString) final String defaultCharsetName = Charset.defaultCharset().name(); final String basePath = URLDecoder.decode(uri.getPath(), defaultCharsetName); - final String query = uri.getQuery(); - final String[] queryKeyValuePairs = query.split("&"); // note: uses "fastpath" for simple regex - String dataSet = null; - Long x = null; - Long y = null; - Long z = null; - Integer width = null; - Integer height = null; - for (final String keyValuePair : queryKeyValuePairs) { - final String[] keyValue = keyValuePair.split("="); - if (keyValue.length == 2) { - final String key = keyValue[0]; - if ("x".equals(key)){ - x = Long.valueOf(keyValue[1]); - } else if ("y".equals(key)) { - y = Long.valueOf(keyValue[1]); - } else if ("z".equals(key)) { - z = Long.valueOf(keyValue[1]); - } else if ("w".equals(key)) { - width = Integer.valueOf(keyValue[1]); - } else if ("h".equals(key)) { - height = Integer.valueOf(keyValue[1]); - } else if ("dataSet".equals(key)) { - dataSet = URLDecoder.decode(keyValue[1], defaultCharsetName); - } - } - } + final QueryKeyValueParameters query = new QueryKeyValueParameters(uri.getQuery(), urlString); - long[] xAndYOffsets = null; - if (x != null) { - if (y != null) { - xAndYOffsets = new long[] { x, y }; - } else { - xAndYOffsets = new long[] { x, 0 }; - } - } else if (y != null) { - xAndYOffsets = new long[] { 0, y }; - } + final Optional rawDataSet = query.getString("dataSet"); + final String dataSet = rawDataSet.isPresent() ? URLDecoder.decode(rawDataSet.get(), defaultCharsetName) : null; + final Long z = query.getLong("z").orElse(null); + Integer width = query.getInt("w").orElse(null); + Integer height = query.getInt("h").orElse(null); + + final Optional x = query.getLong("x"); + final Optional y = query.getLong("y"); + final long[] xAndYOffsets = (x.isPresent() || y.isPresent()) + ? new long[] { x.orElse(0L), y.orElse(0L) } + : null; if ((basePath != null) && (dataSet != null)) { @@ -128,20 +105,13 @@ public ImageProcessor load(final String urlString) height = (int) dimensions[1]; } - switch(dataType) { - case UINT8: - imageProcessor = UNSIGNED_BYTE_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); - break; - case INT16: - imageProcessor = SHORT_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); - break; - case FLOAT32: - imageProcessor = FLOAT_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); - break; - default: - // case INT8: case INT32: case INT64: case FLOAT64: case OBJECT: case UINT16: case UINT32: case UINT64: - throw new IllegalArgumentException("dataType " + dataType + " is not supported"); - } + imageProcessor = switch (dataType) { + case UINT8 -> UNSIGNED_BYTE_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); + case INT16 -> SHORT_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); + case FLOAT32 -> FLOAT_HELPER.load(reader, dataSet, width, height, xAndYOffsets, z); + // case INT8: case INT32: case INT64: case FLOAT64: case OBJECT: case UINT16: case UINT32: case UINT64: + default -> throw new IllegalArgumentException("dataType " + dataType + " is not supported"); + }; } else { throw new IllegalArgumentException( diff --git a/render-app/src/main/java/org/janelia/alignment/transform/DisplacementFieldTransform.java b/render-app/src/main/java/org/janelia/alignment/transform/DisplacementFieldTransform.java new file mode 100644 index 000000000..4390ceb20 --- /dev/null +++ b/render-app/src/main/java/org/janelia/alignment/transform/DisplacementFieldTransform.java @@ -0,0 +1,333 @@ +package org.janelia.alignment.transform; + +import mpicbg.trakem2.transform.CoordinateTransform; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealRandomAccess; +import net.imglib2.RealRandomAccessible; +import net.imglib2.converter.Converters; +import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory; +import net.imglib2.realtransform.AffineTransform2D; +import net.imglib2.realtransform.RealViews; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; +import org.janelia.saalfeldlab.n5.KeyValueAccess; +import org.janelia.saalfeldlab.n5.N5Reader; +import org.janelia.saalfeldlab.n5.googlecloud.GoogleCloudStorageKeyValueAccess; +import org.janelia.saalfeldlab.n5.imglib2.N5Utils; +import org.janelia.n5.precomputed.N5PrecomputedReader; +import org.janelia.n5.precomputed.PrecomputedKeyValueReader; + +import com.google.cloud.storage.Storage; +import com.google.cloud.storage.StorageOptions; +import com.google.gson.GsonBuilder; +import org.janelia.alignment.util.QueryKeyValueParameters; +import org.janelia.saalfeldlab.n5.universe.N5Factory; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.net.URI; +import java.util.Map; +import java.util.Optional; +import java.util.Set; +import java.util.concurrent.ConcurrentHashMap; + + +/** + * Moves each queried location by a displacement vector interpolated from a field on disk. The field is a pull + * map (see {@link #extractAndTransform}), so applying it means inverting the field, done in {@link #applyInPlace} + * by fixed-point iteration. + */ +public class DisplacementFieldTransform + implements CoordinateTransform { + + /** URI (as supplied to {@link #init}) identifying the field on disk and its world-coordinate mapping. */ + private String fieldSourceUri; + private int fieldZIndex; + /** Full-resolution pixels per field pixel; the same in x and y. */ + private double scale; + /** World coordinate that field index 0 maps to, in x and y. */ + private double[] offset; + /** Full-resolution pixels per unit of stored vector. */ + private double vectorScale; + + // ImgLib2 accessor for the displacement field; null until a field source has been loaded. + private RealRandomAccess displacementX; + private RealRandomAccess displacementY; + + /** Whether this instance has already logged a non-converging inversion (see {@link #applyInPlace}). */ + private boolean divergenceLogged = false; + + /** + * Reflection constructor; leaves the instance uninitialized until {@link #init(String)} is called. + */ + public DisplacementFieldTransform() { + this.fieldSourceUri = null; + this.fieldZIndex = -1; + this.scale = DEFAULT_SCALE; + this.offset = new double[] { DEFAULT_OFFSET, DEFAULT_OFFSET }; + this.vectorScale = DEFAULT_VECTOR_SCALE; + + this.displacementX = null; + this.displacementY = null; + } + + /** + * Constructs and immediately loads a transform for the field at the specified source. + * + * @param fieldSourceUri URI locating the field on disk (see class Javadoc for format). + * @param fieldZIndex The z-slice index of the field to use (the field may be 3D, but this transform is 2D) + * @param scale Full-resolution pixels per field pixel, in x and y alike (1 leaves the field at full resolution) + * @param offset World coordinate that field index 0 maps to, in x and y (0 puts the field at the world origin) + * @param vectorScale Full-resolution pixels per unit of stored vector (1 for vectors already in + * full-resolution units, which is what SOFIMA emits) + * + * @throws IllegalArgumentException + * if the field cannot be loaded. + */ + public DisplacementFieldTransform(final String fieldSourceUri, + final int fieldZIndex, + final double scale, + final double[] offset, + final double vectorScale) { + this.init(fieldSourceUri, fieldZIndex, scale, offset, vectorScale); + } + + private void init(final String fieldSourceUri, + final int fieldZIndex, + final double scale, + final double[] offset, + final double vectorScale) { + this.fieldSourceUri = fieldSourceUri; + this.fieldZIndex = fieldZIndex; + this.scale = scale; + this.offset = offset; + this.vectorScale = vectorScale; + + // SOFIMA output as a Neuroglancer precomputed volume, layout [x,y,z,channel], channel 0/1 = X/Y vectors. + final RandomAccessibleInterval fieldRaw = openRawField(fieldSourceUri); + + // x/y out of range is handled by the mirrored extension in extractAndTransform; z is not, so check it here. + if ((fieldZIndex < 0) || (fieldZIndex >= fieldRaw.dimension(2))) { + throw new IllegalArgumentException( + "z " + fieldZIndex + " is outside the z range [0, " + fieldRaw.dimension(2) + + ") of the field at " + fieldSourceUri); + } + + displacementX = extractAndTransform(fieldRaw, 0); + displacementY = extractAndTransform(fieldRaw, 1); + } + + /** + * Cache of raw fields keyed by source URI, so a layer's tiles share one reader and chunk cache instead of each + * re-opening it. Per-instance accessors are still built fresh in {@link #extractAndTransform} since imglib2 + * accessors aren't thread-safe. + */ + private static final Map> RAW_FIELD_CACHE = new ConcurrentHashMap<>(); + + private static RandomAccessibleInterval openRawField(final String fieldSourceUri) { + return RAW_FIELD_CACHE.computeIfAbsent(fieldSourceUri, uri -> { + final N5Reader fieldReader = openPrecomputedReader(uri); + final String scaleKey = fieldReader.list("/")[0]; + return N5Utils.open(fieldReader, scaleKey); + }); + } + + /** + * Opens a Neuroglancer precomputed field (optionally {@code precomputed://}-prefixed) through the N5 API. + * Wired up by hand since {@code n5-universe}'s {@code N5Factory} doesn't know the precomputed format yet. + * {@code gs://} buckets are read anonymously; other schemes go through {@link N5Factory}'s key-value access. + * Exposed so field-preparing clients (e.g. {@code ImportSofimaClient}) use the same path. The dataset lives + * under the first scale key, {@code reader.list("/")[0]}. + */ + public static N5Reader openPrecomputedReader(final String fieldSourceUri) { + String uri = fieldSourceUri; + if (uri.startsWith("precomputed://")) { + uri = uri.substring("precomputed://".length()); + } + + if (uri.startsWith("gs://")) { + // gs:// buckets are read anonymously (the public warp-field bucket needs no credentials). + final Storage storage = StorageOptions.getUnauthenticatedInstance().getService(); + final KeyValueAccess keyValueAccess = new GoogleCloudStorageKeyValueAccess(storage, uri, false); + return new PrecomputedKeyValueReader(keyValueAccess, uri, new GsonBuilder(), true); + } + + // Local filesystem (optionally file://-prefixed): N5PrecomputedReader wires up FileSystemKeyValueAccess + // over the default filesystem. (n5-universe 1.6.0's N5Factory.getKeyValueAccess is package-private.) + final String path = uri.startsWith("file://") ? URI.create(uri).getPath() : uri; + return new N5PrecomputedReader(path, new GsonBuilder(), true); + } + + /** + * The field is a pull map: the vector at a target position points at the source it was pulled from, + * i.e. {@code source = target + vector}. Render's transform lists run source to target, so the vectors are + * negated here and scaled by {@code vectorScale} to full resolution. + */ + private RealRandomAccess extractAndTransform(final RandomAccessibleInterval rawField, + final int xory) { + // Replace NaNs before interpolating, so they don't leak into neighboring pixels. + final RandomAccessibleInterval cleaned = Converters.convertRAI( + rawField, + (i, o) -> o.set(Float.isNaN(i.getRealFloat()) ? 0 : i.getRealFloat()), + new FloatType()); + + // Slice the [x,y,z,channel] dataset: choose the vector component (channel, dim=3) and then the + // z-slice (dim=2). Slicing the higher dimension (channel) first keeps the z index valid at dim=2. + final RandomAccessibleInterval slice = Views.hyperSlice( + Views.hyperSlice(cleaned, 3, xory), 2, this.fieldZIndex); + + // Place the slice in world coordinates: field index 0 lands on offset, one field pixel spans scale + // full-resolution pixels, so a query at p reads the field at (p - offset) / scale. + final AffineTransform2D fieldToWorld = new AffineTransform2D(); + fieldToWorld.set(this.scale, 0, this.offset[0], + 0, this.scale, this.offset[1]); + final RealRandomAccessible scaledAndInterpolated = RealViews.affine( + Views.interpolate(Views.extendMirrorDouble(slice), new NLinearInterpolatorFactory<>()), + fieldToWorld); + + // Negate and scale in one pass, applied after interpolation. Negating just flips the vectors; the actual + // inversion (evaluating at the target rather than the source) happens in applyInPlace. + final double pullToPushScale = -this.vectorScale; + return Converters.convert( + scaledAndInterpolated, + (i, o) -> o.set((float) (i.getRealFloat() * pullToPushScale)), + new FloatType()).realRandomAccess(); + } + + @Override + public double[] apply(final double[] location) { + final double[] out = location.clone(); + applyInPlace(out); + return out; + } + + @Override + public void applyInPlace(final double[] location) { + + if (displacementX == null || displacementY == null) { + throw new IllegalStateException( + "displacement field has not been loaded; call init(String) before applying this transform"); + } + + // The (negated) field vector belongs to the target location, not to the queried source location, so the + // target solves t = source + d(t). Iterate t <- source + d(t) from t = source; this converges (linearly, + // at the rate of the field's Jacobian norm) for any field whose Jacobian norm stays below one, which is + // the same condition under which the field is invertible at all. + final double[] target = { location[0], location[1] }; + final double[] vector = new double[2]; + double step = Double.POSITIVE_INFINITY; + for (int i = 0; (i < MAX_INVERSION_ITERATIONS) && (step > INVERSION_TOLERANCE); i++) { + lookUpVector(target, vector); + final double x = location[0] + vector[0]; + final double y = location[1] + vector[1]; + step = Math.max(Math.abs(x - target[0]), Math.abs(y - target[1])); + target[0] = x; + target[1] = y; + } + + // A non-invertible field may still exceed the tolerance after the cap; use the last estimate rather than + // failing the render, logged once per instance to avoid a line per pixel. + if ((step > INVERSION_TOLERANCE) && (! divergenceLogged)) { + divergenceLogged = true; + LOG.warn("applyInPlace: inversion did not converge within {} iterations at ({}, {}) for field {}; " + + "residual step is {} px and the last estimate is used. Further occurrences for this " + + "transform instance are not logged.", + MAX_INVERSION_ITERATIONS, location[0], location[1], toDataString(), step); + } + + location[0] = target[0]; + location[1] = target[1]; + } + + /** + * Looks up the interpolated field vector (already negated and scaled to full resolution) at a world location. + * Package private so that tests can check the field placement on its own, separately from the inversion. + */ + void lookUpVector(final double[] location, final double[] vector) { + vector[0] = displacementX.setPositionAndGet(location).getRealDouble(); + vector[1] = displacementY.setPositionAndGet(location).getRealDouble(); + } + + /** + * Parses the data string (field source URI plus {@code ?key=value} params, e.g. + * {@code file:///path/to/field.n5?z=5&scale=40.0&offset=-5318.0,-783.0}) and loads the field. Only + * {@code z} is required; the rest default to the identity placement. Unknown parameters are rejected so a + * misspelled one can't silently fall back to its default. + * + * @throws IllegalArgumentException + * if the data string cannot be parsed or the field cannot be loaded. + */ + @Override + public void init(final String data) throws IllegalArgumentException { + + final String trimmed = data.trim(); + final int queryStart = trimmed.indexOf('?'); + if (queryStart < 0) { + throw new IllegalArgumentException( + "transform data must be a field source URI followed by '?z=' and optionally " + + "'&scale=&offset=,&vectorScale=', " + + "but was '" + data + "'"); + } + + final String parsedSourceUri = trimmed.substring(0, queryStart); + final QueryKeyValueParameters params = new QueryKeyValueParameters(trimmed.substring(queryStart + 1), data); + params.validateKeys(VALID_PARAMETERS); + + final Optional parsedOffset = params.getDoubleArray("offset"); + if (parsedOffset.isPresent() && (parsedOffset.get().length != 2)) { + throw new IllegalArgumentException( + "parameter 'offset' must be two comma separated numbers, but was '" + + params.getString("offset").orElseThrow() + "' in transform data '" + data + "'"); + } + final double[] offset = parsedOffset.orElse(new double[] { DEFAULT_OFFSET, DEFAULT_OFFSET }); + + init(parsedSourceUri, + params.getInt("z").orElseThrow(() -> new IllegalArgumentException("missing required parameter 'z' in transform data '" + data + "'")), + params.getDouble("scale").orElse(DEFAULT_SCALE), + offset, + params.getDouble("vectorScale").orElse(DEFAULT_VECTOR_SCALE)); + } + + @Override + public String toXML(final String indent) { + return indent + ""; + } + + @Override + public String toDataString() { + // Writes every parameter, even defaults, so a persisted string keeps its meaning if a default ever changes. + return fieldSourceUri + + "?z=" + fieldZIndex + + "&scale=" + scale + + "&offset=" + offset[0] + "," + offset[1] + + "&vectorScale=" + vectorScale; + } + + @Override + public CoordinateTransform copy() { + // Re-loads the field so the copy has independent accessors (imglib2 accessors are not thread-safe). + return new DisplacementFieldTransform(fieldSourceUri, fieldZIndex, scale, offset.clone(), vectorScale); + } + + @Override + public String toString() { + return "{ \"fieldSourceUri\": \"" + fieldSourceUri + + "\", \"fieldZIndex\": " + fieldZIndex + + ", \"scale\": " + scale + + ", \"offset\": [" + offset[0] + ", " + offset[1] + "]" + + ", \"vectorScale\": " + vectorScale + " }"; + } + + private static final Logger LOG = LoggerFactory.getLogger(DisplacementFieldTransform.class); + + private static final double DEFAULT_SCALE = 1.0; + private static final double DEFAULT_OFFSET = 0.0; + private static final double DEFAULT_VECTOR_SCALE = 1.0; + + /** Full-resolution pixels of movement below which the inversion in {@link #applyInPlace} is considered done. */ + private static final double INVERSION_TOLERANCE = 1e-4; + private static final int MAX_INVERSION_ITERATIONS = 20; + + private static final Set VALID_PARAMETERS = Set.of("z", "scale", "offset", "vectorScale"); +} diff --git a/render-app/src/main/java/org/janelia/alignment/util/QueryKeyValueParameters.java b/render-app/src/main/java/org/janelia/alignment/util/QueryKeyValueParameters.java new file mode 100644 index 000000000..821871f67 --- /dev/null +++ b/render-app/src/main/java/org/janelia/alignment/util/QueryKeyValueParameters.java @@ -0,0 +1,173 @@ +package org.janelia.alignment.util; + +import java.util.HashMap; +import java.util.Map; +import java.util.Optional; +import java.util.Set; + +/** + * Parses a {@code key=value&key2=value2} style query string into typed accessors. + * Used both for URI query strings and for the {@code ?key=value&...} suffix of a + * hand-encoded data string (e.g. a {@code CoordinateTransform} data string). + *

+ * Every accessor returns an empty {@link Optional} for an absent key; callers decide how to handle + * that (a required parameter throws via {@code orElseThrow}, an optional one falls back via + * {@code orElse}). + */ +public class QueryKeyValueParameters { + + private final Map params = new HashMap<>(); + private final String source; + + /** + * @param query the query string to parse (everything after the {@code ?}, without the {@code ?} itself); + * may be {@code null} or empty, in which case no parameters are present. + * @param source the full original string, included in error messages so they can point back to it. + */ + public QueryKeyValueParameters(final String query, final String source) { + this.source = source; + if (query != null) { + for (final String pair : query.split("&")) { + if (pair.isEmpty()) { + continue; + } + final int eq = pair.indexOf('='); + if (eq < 0) { + throw new IllegalArgumentException( + "invalid query parameter '" + pair + "' in '" + source + "'"); + } + params.put(pair.substring(0, eq), pair.substring(eq + 1)); + } + } + } + + public Optional getString(final String key) { + return Optional.ofNullable(params.get(key)); + } + + /** @throws IllegalArgumentException if present but not a valid integer. */ + public Optional getInt(final String key) { + return getString(key).map(value -> parseInt(key, value)); + } + + /** @throws IllegalArgumentException if present but not a valid long. */ + public Optional getLong(final String key) { + return getString(key).map(value -> parseLong(key, value)); + } + + /** @throws IllegalArgumentException if present but not a valid float. */ + public Optional getFloat(final String key) { + return getString(key).map(value -> parseFloat(key, value)); + } + + /** @throws IllegalArgumentException if present but not a valid double. */ + public Optional getDouble(final String key) { + return getString(key).map(value -> parseDouble(key, value)); + } + + /** @throws IllegalArgumentException if present but contains an invalid integer. */ + public Optional getIntArray(final String key) { + return getString(key).map(value -> parseIntArray(key, value)); + } + + /** @throws IllegalArgumentException if present but contains an invalid long. */ + public Optional getLongArray(final String key) { + return getString(key).map(value -> parseLongArray(key, value)); + } + + /** @throws IllegalArgumentException if present but contains an invalid float. */ + public Optional getFloatArray(final String key) { + return getString(key).map(value -> parseFloatArray(key, value)); + } + + /** @throws IllegalArgumentException if present but contains an invalid double. */ + public Optional getDoubleArray(final String key) { + return getString(key).map(value -> parseDoubleArray(key, value)); + } + + /** + * @throws IllegalArgumentException if any parsed key is not in {@code validKeys}; helps catch typos that + * would otherwise silently fall back to a default value. + */ + public void validateKeys(final Set validKeys) { + for (final String key : params.keySet()) { + if (! validKeys.contains(key)) { + throw new IllegalArgumentException( + "unknown query parameter '" + key + "' in '" + source + + "'; supported parameters are " + validKeys); + } + } + } + + private int parseInt(final String key, final String value) { + try { + return Integer.parseInt(value); + } catch (final NumberFormatException e) { + throw new IllegalArgumentException( + "invalid integer value '" + value + "' for parameter '" + key + "' in '" + source + "'", e); + } + } + + private long parseLong(final String key, final String value) { + try { + return Long.parseLong(value); + } catch (final NumberFormatException e) { + throw new IllegalArgumentException( + "invalid long value '" + value + "' for parameter '" + key + "' in '" + source + "'", e); + } + } + + private float parseFloat(final String key, final String value) { + try { + return Float.parseFloat(value); + } catch (final NumberFormatException e) { + throw new IllegalArgumentException( + "invalid float value '" + value + "' for parameter '" + key + "' in '" + source + "'", e); + } + } + + private double parseDouble(final String key, final String value) { + try { + return Double.parseDouble(value); + } catch (final NumberFormatException e) { + throw new IllegalArgumentException( + "invalid double value '" + value + "' for parameter '" + key + "' in '" + source + "'", e); + } + } + + private int[] parseIntArray(final String key, final String value) { + final String[] parts = value.split(","); + final int[] result = new int[parts.length]; + for (int i = 0; i < parts.length; i++) { + result[i] = parseInt(key, parts[i]); + } + return result; + } + + private long[] parseLongArray(final String key, final String value) { + final String[] parts = value.split(","); + final long[] result = new long[parts.length]; + for (int i = 0; i < parts.length; i++) { + result[i] = parseLong(key, parts[i]); + } + return result; + } + + private float[] parseFloatArray(final String key, final String value) { + final String[] parts = value.split(","); + final float[] result = new float[parts.length]; + for (int i = 0; i < parts.length; i++) { + result[i] = parseFloat(key, parts[i]); + } + return result; + } + + private double[] parseDoubleArray(final String key, final String value) { + final String[] parts = value.split(","); + final double[] result = new double[parts.length]; + for (int i = 0; i < parts.length; i++) { + result[i] = parseDouble(key, parts[i]); + } + return result; + } +} diff --git a/render-app/src/test/java/org/janelia/alignment/transform/DisplacementFieldTransformTest.java b/render-app/src/test/java/org/janelia/alignment/transform/DisplacementFieldTransformTest.java new file mode 100644 index 000000000..b6246b988 --- /dev/null +++ b/render-app/src/test/java/org/janelia/alignment/transform/DisplacementFieldTransformTest.java @@ -0,0 +1,185 @@ +package org.janelia.alignment.transform; + +import java.nio.file.Path; +import java.util.Arrays; + +import org.janelia.saalfeldlab.n5.DataType; +import org.janelia.n5.precomputed.PrecomputedTestVolumes; +import org.junit.Assert; +import org.junit.Rule; +import org.junit.Test; +import org.junit.rules.TemporaryFolder; + +import mpicbg.trakem2.transform.CoordinateTransform; + +/** + * Tests the {@link DisplacementFieldTransform} class. + */ +public class DisplacementFieldTransformTest { + + private static final String SAMPLE_URI = + "file:///tmp/does-not-exist.n5?z=5&scale=8.0&offset=100.0,-50.0&vectorScale=2.0"; + + @Test + public void testDataStringRoundTrip() { + // init should fail to open the (nonexistent) field, but only after parsing the data string, + // so a real field is not needed to exercise the parse-and-serialize contract. + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + try { + transform.init(SAMPLE_URI); + Assert.fail("expected init to fail loading a nonexistent field"); + } catch (final RuntimeException e) { + Assert.assertEquals("data string should round-trip even when loading fails", + SAMPLE_URI, transform.toDataString()); + } + } + + @Test + public void testApplyBeforeInitFails() { + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + try { + transform.applyInPlace(new double[] {0.0, 0.0}); + Assert.fail("expected applyInPlace to fail before the field is loaded"); + } catch (final IllegalStateException e) { + Assert.assertTrue("exception should mention initialization", + e.getMessage().contains("init")); + } + } + + @Rule + public TemporaryFolder tempFolder = new TemporaryFolder(); + + /** + * End-to-end read of a real Neuroglancer-precomputed field through n5-ng-precomputed: guards the + * [x,y,z,channel] slicing, both scalings, and (since it actually touches the reader) the n5 version + * alignment between this reactor and the locally installed n5-ng-precomputed build. + */ + @Test + public void testAppliesPrecomputedField() throws Exception { + + // value at (x,y,z) is x+y+10*z for the X component (channel 0) and 100 more for Y (channel 1) + final Path fieldDir = tempFolder.newFolder("field").toPath(); + PrecomputedTestVolumes.writeRawVolume(fieldDir, + DataType.FLOAT32, + 2, + new long[] {4, 3, 2}, + new int[] {4, 3, 2}, + new long[] {0, 0, 0}, + (x, y, z, c) -> x + y + 10 * z + 100 * c); + + // defaults only: scale 1 and offset 0 mean (1,1) reads field (1,1) of z-slice 1, and vector scale 1 means + // the stored vectors are used as-is apart from the pull-to-push negation + assertDisplacement(fieldDir, "z=1", new double[] {1.0, 1.0}, -12.0, -112.0); + + // scale 2 halves the query position (so (2,2) reads field (1,1)) and vectorScale 4 quadruples the vectors + assertDisplacement(fieldDir, "z=1&scale=2.0&vectorScale=4.0", + new double[] {2.0, 2.0}, -12.0 * 4, -112.0 * 4); + + // offset shifts the query position, so (2,2) again reads field (1,1) + assertDisplacement(fieldDir, "z=1&offset=1.0,1.0", + new double[] {2.0, 2.0}, -12.0, -112.0); + + // a malformed offset is rejected rather than silently read as one number + final DisplacementFieldTransform malformed = new DisplacementFieldTransform(); + try { + malformed.init(fieldDir + "?z=1&offset=1.0"); + Assert.fail("expected init to reject a one-component offset"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue("exception should mention the offset, but was: " + e.getMessage(), + e.getMessage().contains("offset")); + } + + // x and y beyond the field are answered from the mirrored extension rather than failing; just past the last + // sample (the field is 4 wide, so x=3) the double-mirrored extension repeats that boundary value + assertDisplacement(fieldDir, "z=1", new double[] {4.0, 1.0}, -14.0, -114.0); + + // z is guarded instead, since an out-of-range slice would read outside the cached image + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + try { + transform.init(fieldDir + "?z=2"); + Assert.fail("expected init to reject a z index outside the field"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue("exception should mention the z range, but was: " + e.getMessage(), + e.getMessage().contains("z range")); + } + } + + private static void assertDisplacement(final Path fieldDir, + final String queryString, + final double[] location, + final double expectedDx, + final double expectedDy) { + + final String data = fieldDir + "?" + queryString; + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + transform.init(data); + + // the raw lookup, not apply(): apply() inverts the field, which would move the query off the position + // whose placement is under test here (and this steep test field is not invertible anyway) + final double[] vector = new double[2]; + transform.lookUpVector(location, vector); + Assert.assertEquals("wrong x displacement for " + data, expectedDx, vector[0], 0.0001); + Assert.assertEquals("wrong y displacement for " + data, expectedDy, vector[1], 0.0001); + } + + /** + * The field is a pull map, so applying it means solving {@code t = p + d(t)} rather than evaluating {@code d} + * at {@code p}. Uses a field with displacement {@code -0.2*(x,y)}, whose exact fixed point is {@code p/1.2} — + * clearly apart from the first-order answer {@code 0.8*p}. + */ + @Test + public void testInvertsFieldByIteration() throws Exception { + + final Path fieldDir = tempFolder.newFolder("rampField").toPath(); + PrecomputedTestVolumes.writeRawVolume(fieldDir, + DataType.FLOAT32, + 2, + new long[] {64, 48, 1}, + new int[] {64, 48, 1}, + new long[] {0, 0, 0}, + (x, y, z, c) -> (c == 0) ? x : y); + + // vectorScale carries the 0.2 because the test volume can only hold whole numbers + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + transform.init(fieldDir + "?z=0&vectorScale=0.2"); + + final double[] target = transform.apply(new double[] {24.0, 12.0}); + Assert.assertEquals("x should solve t = 24 - 0.2 * t", 20.0, target[0], 0.001); + Assert.assertEquals("y should solve t = 12 - 0.2 * t", 10.0, target[1], 0.001); + + // the defining property, independent of the analytic solution above + final double[] vector = new double[2]; + transform.lookUpVector(target, vector); + Assert.assertEquals("x residual", 0.0, 24.0 + vector[0] - target[0], 0.001); + Assert.assertEquals("y residual", 0.0, 12.0 + vector[1] - target[1], 0.001); + + // the same field at vectorScale 2 is not invertible (Jacobian norm 2), so the iteration hits its cap: + // that must warn and return the last estimate rather than throw or hand back a non-number + final DisplacementFieldTransform steep = new DisplacementFieldTransform(); + steep.init(fieldDir + "?z=0&vectorScale=2.0"); + final double[] estimate = steep.apply(new double[] {24.0, 12.0}); + Assert.assertTrue("a non-converging inversion should still return numbers, but was " + + Arrays.toString(estimate), + Double.isFinite(estimate[0]) && Double.isFinite(estimate[1])); + } + + @Test + public void testMisspelledParameterFails() { + // since everything but z is optional, a typo would otherwise silently apply the default + final DisplacementFieldTransform transform = new DisplacementFieldTransform(); + try { + transform.init("file:///tmp/does-not-exist.n5?z=0&scalex=40.0"); + Assert.fail("expected init to reject the misspelled parameter"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue("exception should name the offending parameter", + e.getMessage().contains("scalex")); + } + } + + @Test + public void testImplementsCoordinateTransform() { + // guards the reflective LeafTransformSpec.newInstance() contract (no-arg constructor + interface) + final CoordinateTransform transform = new DisplacementFieldTransform(); + Assert.assertNotNull("no-arg constructed instance should exist", transform); + } +} diff --git a/render-app/src/test/java/org/janelia/alignment/util/QueryKeyValueParametersTest.java b/render-app/src/test/java/org/janelia/alignment/util/QueryKeyValueParametersTest.java new file mode 100644 index 000000000..0994dc589 --- /dev/null +++ b/render-app/src/test/java/org/janelia/alignment/util/QueryKeyValueParametersTest.java @@ -0,0 +1,97 @@ +package org.janelia.alignment.util; + +import java.util.Optional; +import java.util.Set; + +import org.junit.Assert; +import org.junit.Test; + +/** + * Tests the {@link QueryKeyValueParameters} class. + */ +public class QueryKeyValueParametersTest { + + private static final String SOURCE = "test-source"; + + @Test + public void testParsesEveryType() { + final QueryKeyValueParameters params = new QueryKeyValueParameters( + "s=hello&i=42&l=9999999999&f=1.5&d=2.5&ia=1,2,3&la=1,2,3&fa=1.5,2.5&da=1.5,2.5", + SOURCE); + + Assert.assertEquals("hello", params.getString("s").orElseThrow()); + Assert.assertEquals(Integer.valueOf(42), params.getInt("i").orElseThrow()); + Assert.assertEquals(Long.valueOf(9999999999L), params.getLong("l").orElseThrow()); // exceeds int range + Assert.assertEquals(1.5f, params.getFloat("f").orElseThrow(), 0.0001f); + Assert.assertEquals(2.5, params.getDouble("d").orElseThrow(), 0.0001); + Assert.assertArrayEquals(new int[] {1, 2, 3}, params.getIntArray("ia").orElseThrow()); + Assert.assertArrayEquals(new long[] {1, 2, 3}, params.getLongArray("la").orElseThrow()); + Assert.assertArrayEquals(new float[] {1.5f, 2.5f}, params.getFloatArray("fa").orElseThrow(), 0.0001f); + Assert.assertArrayEquals(new double[] {1.5, 2.5}, params.getDoubleArray("da").orElseThrow(), 0.0001); + } + + @Test + public void testAbsentKeyReturnsEmptyForEveryType() { + final QueryKeyValueParameters params = new QueryKeyValueParameters("present=1", SOURCE); + + Assert.assertEquals(Optional.empty(), params.getString("missing")); + Assert.assertEquals(Optional.empty(), params.getInt("missing")); + Assert.assertEquals(Optional.empty(), params.getLong("missing")); + Assert.assertEquals(Optional.empty(), params.getFloat("missing")); + Assert.assertEquals(Optional.empty(), params.getDouble("missing")); + Assert.assertEquals(Optional.empty(), params.getIntArray("missing")); + Assert.assertEquals(Optional.empty(), params.getLongArray("missing")); + Assert.assertEquals(Optional.empty(), params.getFloatArray("missing")); + Assert.assertEquals(Optional.empty(), params.getDoubleArray("missing")); + } + + @Test + public void testMalformedPairFails() { + try { + new QueryKeyValueParameters("noEqualsSign", SOURCE); + Assert.fail("expected a pair without '=' to fail"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue(e.getMessage().contains("noEqualsSign")); + } + } + + @Test + public void testInvalidNumberFailsForEveryScalarType() { + final QueryKeyValueParameters params = new QueryKeyValueParameters("v=notANumber", SOURCE); + assertInvalid(() -> params.getInt("v"), "integer"); + assertInvalid(() -> params.getLong("v"), "long"); + assertInvalid(() -> params.getFloat("v"), "float"); + assertInvalid(() -> params.getDouble("v"), "double"); + } + + @Test + public void testInvalidNumberFailsForEveryArrayType() { + final QueryKeyValueParameters params = new QueryKeyValueParameters("v=1,notANumber", SOURCE); + assertInvalid(() -> params.getIntArray("v"), "integer"); + assertInvalid(() -> params.getLongArray("v"), "long"); + assertInvalid(() -> params.getFloatArray("v"), "float"); + assertInvalid(() -> params.getDoubleArray("v"), "double"); + } + + private static void assertInvalid(final Runnable action, final String expectedTypeName) { + try { + action.run(); + Assert.fail("expected invalid " + expectedTypeName + " value to fail"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue("exception should mention " + expectedTypeName + ", but was: " + e.getMessage(), + e.getMessage().contains(expectedTypeName)); + } + } + + @Test + public void testValidateKeysRejectsUnknownKey() { + final QueryKeyValueParameters params = new QueryKeyValueParameters("a=1&b=2", SOURCE); + try { + params.validateKeys(Set.of("a")); + Assert.fail("expected unknown key 'b' to fail validation"); + } catch (final IllegalArgumentException e) { + Assert.assertTrue(e.getMessage().contains("b")); + } + params.validateKeys(Set.of("a", "b")); // should not throw once every key is known + } +} diff --git a/render-ws-java-client/src/main/java/org/janelia/render/client/multisem/ImportSofimaClient.java b/render-ws-java-client/src/main/java/org/janelia/render/client/multisem/ImportSofimaClient.java new file mode 100644 index 000000000..bfedbd254 --- /dev/null +++ b/render-ws-java-client/src/main/java/org/janelia/render/client/multisem/ImportSofimaClient.java @@ -0,0 +1,206 @@ +package org.janelia.render.client.multisem; + +import com.beust.jcommander.Parameter; +import com.beust.jcommander.ParametersDelegate; + +import java.io.IOException; +import java.util.Arrays; +import java.util.List; +import java.util.concurrent.ForkJoinPool; +import java.util.concurrent.ForkJoinTask; + +import org.janelia.alignment.spec.Bounds; +import org.janelia.alignment.spec.LeafTransformSpec; +import org.janelia.alignment.spec.ResolvedTileSpecCollection; +import org.janelia.alignment.spec.stack.StackMetaData; +import org.janelia.alignment.transform.DisplacementFieldTransform; +import org.janelia.render.client.ClientRunner; +import org.janelia.render.client.RenderDataClient; +import org.janelia.render.client.parameter.CommandLineParameters; +import org.janelia.render.client.parameter.RenderWebServiceParameters; +import org.janelia.render.client.parameter.ZRangeParameters; +import org.janelia.saalfeldlab.n5.N5Reader; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +/** + * Adds a {@link DisplacementFieldTransform} to every tile spec of a stack, layer by layer. + *

+ * The (dense) displacement field is currently expected to be SOFIMA output for multi-SEM acquisitions, stored as a + * Neuroglancer precomputed volume with a 4D {@code [x,y,z,channel]} layout and one z-slice per stack layer (opened + * through the same {@link DisplacementFieldTransform#openPrecomputedReader} path the transform itself uses). For + * each layer, this client + *

    + *
  • computes {@code fieldZIndex} from the layer z and the stack's {@code minZ},
  • + *
  • appends a {@link DisplacementFieldTransform} with the resulting data string to each tile spec, and
  • + *
  • saves the modified tile specs to the target stack.
  • + *
+ * Field index {@code (0,0,0)} is taken to sit on the minimum corner {@code (minX,minY,minZ)} of the source stack + * bounds, since the field is computed on an export of that stack and an export re-origins the data at {@code (0,0,0)}. + * The x and y parts of that corner go into the data string as the transform's offset; the z part turns the layer z + * into the field's z-slice index. The data string's scale is either the {@code --scale} given on the command line or, + * if that is omitted, the stack bounds divided by the field dimensions and rounded to a whole number. Only the vector + * scale is left at its default, which suits SOFIMA output with vectors already in full-resolution units. + * Layers are processed in order, but the tiles within a layer are processed by {@code --numThreads} threads, which is + * what parallelizes the field chunk reads that deriving the bounding boxes triggers. + * If {@code --targetStack} is given, the modified tile specs are written there (the stack is derived from the source + * if it does not yet exist); otherwise they are written back into the source stack. If {@code --completeTargetStack} + * is set, the target stack is completed once all layers have been saved. + * + * @author Michael Innerberger + */ +public class ImportSofimaClient { + + private final Parameters params; + private final RenderDataClient renderClient; + + public static class Parameters extends CommandLineParameters { + @ParametersDelegate + private final RenderWebServiceParameters renderParams = new RenderWebServiceParameters(); + @ParametersDelegate + private final ZRangeParameters zRangeParams = new ZRangeParameters(); + @Parameter(names = "--stack", description = "Source stack to which the displacement field is added", required = true) + private String stack; + @Parameter(names = "--targetStack", description = "Stack to save modified tile specs to", required = true) + private String targetStack; + @Parameter(names = "--sofimaFieldUri", description = "URI of the SOFIMA displacement field N5 container", required = true) + private String sofimaFieldUri; + @Parameter(names = "--scale", description = "Full-resolution pixels per field pixel, i.e. the factor by which the field is downsampled in x and y (e.g. 40); derived from the stack bounds and the field dimensions if omitted") + private Double scale; + @Parameter(names = "--completeTargetStack", description = "Complete the target stack after all layers have been saved") + private boolean completeTargetStack = false; + @Parameter(names = "--numThreads", description = "Number of tiles within a layer to process concurrently (default: 1)") + private int numThreads = 1; + } + + public static void main(final String[] args) { + final ClientRunner clientRunner = new ClientRunner(args) { + @Override + public void runClient(final String[] args) throws Exception { + final Parameters parameters = new Parameters(); + parameters.parse(args); + LOG.info("runClient: entry, parameters={}", parameters); + + final ImportSofimaClient client = new ImportSofimaClient(parameters); + client.addDisplacementField(); + } + }; + clientRunner.run(); + } + + public ImportSofimaClient(final Parameters parameters) { + this.params = parameters; + this.renderClient = new RenderDataClient(parameters.renderParams.baseDataUrl, + parameters.renderParams.owner, + parameters.renderParams.project); + } + + public void addDisplacementField() throws Exception { + + final StackMetaData sourceStackMetaData = renderClient.getStackMetaData(params.stack); + final Bounds stackBounds = sourceStackMetaData.getStats().getStackBounds(); + + // The field is computed on an export of this stack, and an export puts its own origin at (0,0,0) and merely + // notes the world offset in its metadata - which the field producer does not read. So field index (0,0,0) + // sits on the minimum corner of the stack bounding box, in z just as much as in x and y. + final double[] offset = { stackBounds.getMinX(), stackBounds.getMinY(), stackBounds.getMinZ() }; + + // Open the field up front so that a bad URI fails before any stack is touched, and work out the scale + final double scale; + try (final N5Reader fieldReader = DisplacementFieldTransform.openPrecomputedReader(params.sofimaFieldUri)) { + // The precomputed dataset lives under the first scale key (see DisplacementFieldTransform); the + // layout is [x,y,z,channel], so dim 0 is X and dim 1 is Y. + final String scaleKey = fieldReader.list("/")[0]; + final long[] fieldDimensions = fieldReader.getDatasetAttributes(scaleKey).getDimensions(); + + // The field does not cover the stack bounds exactly, rounding recovers the intended + // factor; the leftover strip is handled by the transform's mirrored extension. + final double xScale = Math.round(stackBounds.getDeltaX() / fieldDimensions[0]); + final double yScale = Math.round(stackBounds.getDeltaY() / fieldDimensions[1]); + if ((params.scale == null) && (xScale != yScale)) { + // The transform downsamples x and y by the same factor, so a field that does not is not supported. + throw new IllegalArgumentException( + "derived x and y scales differ (" + xScale + " vs " + yScale + "); pass --scale explicitly"); + } + scale = (params.scale != null) ? params.scale : xScale; + + LOG.info("addDisplacementField: stack bounds are {}, field {} has dimensions {}, scale is {}, offset is {}", + stackBounds, scaleKey, Arrays.toString(fieldDimensions), scale, Arrays.toString(offset)); + } catch (final Exception e) { + throw new IllegalArgumentException("Failed to process SOFIMA field at " + params.sofimaFieldUri, e); + } + + // Set up the target stack + final String targetStack = params.targetStack; + if (! targetStack.equals(params.stack)) { + renderClient.setupDerivedStack(sourceStackMetaData, targetStack); + } else { + renderClient.ensureStackIsInLoadingState(targetStack, sourceStackMetaData); + } + + // Get and process all z values + final List zValues = renderClient.getStackZValues(params.stack, + params.zRangeParams.minZ, + params.zRangeParams.maxZ); + LOG.info("addDisplacementField: processing {} layers with {} threads", zValues.size(), params.numThreads); + + // One pool, reused for the tiles of each layer in turn (see addFieldToLayer for why tiles and not layers) + try (final ForkJoinPool pool = new ForkJoinPool(params.numThreads)) { + for (final Double z : zValues) { + // Derive the slice from z itself rather than from the running layer index, so that a gap in the + // stack's z values does not shift every later layer onto the wrong slice. + final long fieldZIndex = Math.round(z - offset[2]); + addFieldToLayer(z, buildDataString(fieldZIndex, scale, offset), targetStack, pool); + } + } + + // Complete the target stack + if (params.completeTargetStack) { + LOG.info("addDisplacementField: completing stack {}", targetStack); + renderClient.setStackState(targetStack, StackMetaData.StackState.COMPLETE); + } + + LOG.info("addDisplacementField: exit"); + } + + /** + * Adds the transform to every tile spec of one layer and saves them. Deriving a bounding box evaluates the + * transform over the tile's footprint, so this is where the (blocking) field chunk reads happen. Tiles are the + * natural axis to parallelize: they cover disjoint parts of the field, so concurrent tiles load disjoint chunks + * and the round-trip latency of those reads overlaps. + */ + private void addFieldToLayer(final Double z, + final String dataString, + final String targetStack, + final ForkJoinPool pool) + throws IOException { + + final ResolvedTileSpecCollection tileSpecs = renderClient.getResolvedTiles(params.stack, z); + + // Run the parallel stream inside the given pool. Balances the load internally + pool.invoke(ForkJoinTask.adapt(() -> tileSpecs.getTileSpecs().parallelStream().forEach(tileSpec -> { + final LeafTransformSpec transformSpec = new LeafTransformSpec(DisplacementFieldTransform.class.getName(), dataString); + tileSpec.addTransformSpecs(List.of(transformSpec)); + tileSpec.deriveBoundingBox(tileSpec.getMeshCellSize(), true); + }))); + + renderClient.saveResolvedTiles(tileSpecs, targetStack, z); + LOG.info("addFieldToLayer: saved {} tile specs for z {}", tileSpecs.getTileCount(), z); + } + + /** + * Compiles the {@link DisplacementFieldTransform} data string for one layer. The format must match what + * {@link DisplacementFieldTransform#init(String)} parses. Only the vector scale is omitted, so that the + * transform's default of 1 applies (SOFIMA vectors are already in full-resolution pixels). + */ + private String buildDataString(final long fieldZIndex, + final double scale, + final double[] offset) { + return params.sofimaFieldUri + + "?z=" + fieldZIndex + + "&scale=" + scale + + "&offset=" + offset[0] + "," + offset[1]; + } + + private static final Logger LOG = LoggerFactory.getLogger(ImportSofimaClient.class); +}