From c7044c893bf1e7914e62064ffad4be8dd73f945e Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sat, 11 Jul 2026 13:47:37 -0400 Subject: [PATCH 01/15] Create first version of updated inpainter --- .../spark/multisem/Wafer6061Inpainter.java | 889 ++++++++++++------ 1 file changed, 609 insertions(+), 280 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index b05bc3090..f490d66cf 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -2,15 +2,33 @@ import com.beust.jcommander.Parameter; + +import java.io.IOException; +import java.io.Serializable; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Iterator; +import java.util.LinkedHashMap; +import java.util.List; +import java.util.Map; +import java.util.function.Consumer; + +import bdv.export.Downsample; import net.imglib2.Cursor; -import net.imglib2.Interval; +import net.imglib2.IterableInterval; +import net.imglib2.KDTree; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealPoint; import net.imglib2.img.Img; import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.AccessFlags; +import net.imglib2.img.cell.CellGrid; +import net.imglib2.neighborsearch.KNearestNeighborSearchOnKDTree; import net.imglib2.type.numeric.integer.UnsignedByteType; -import net.imglib2.util.Intervals; +import net.imglib2.type.numeric.real.DoubleType; import net.imglib2.view.Views; + import org.apache.spark.SparkConf; import org.apache.spark.api.java.JavaRDD; import org.apache.spark.api.java.JavaSparkContext; @@ -19,6 +37,7 @@ import org.janelia.render.client.ClientRunner; import org.janelia.render.client.parameter.CommandLineParameters; import org.janelia.render.client.spark.LogUtilities; +import org.janelia.saalfeldlab.n5.DataBlock; import org.janelia.saalfeldlab.n5.DatasetAttributes; import org.janelia.saalfeldlab.n5.N5Reader; import org.janelia.saalfeldlab.n5.N5Writer; @@ -27,20 +46,23 @@ import org.slf4j.Logger; import org.slf4j.LoggerFactory; -import java.io.Serializable; -import java.util.Arrays; -import java.util.List; -import java.util.Objects; - /** - * Class for inpainting small gaps between tiles in the wafer 60/61 dataset. + * Spark client for filling holes in the tissue of wafer 60/61 N5 volumes. *

- * The regions to inpaint are determined by looking at mask pixels: if an unmasked pixel is encountered, - * the algorithm looks at pairs of pixels `stepSize` away in the x and y directions to find non-masked pixels. - * If such a pair is found, the pixel is most likely in a narrow gap between tiles and needs inpainting. The - * inpainting is done by averaging the image data from the adjacent pixels in the z direction. If only one of the - * image values in z is available, that value is used. If neither is available, the pixel is set to 0. + * Inputs are (a) an N5 tissue volume (only non-empty blocks stored), (b) an N5 mask marking where image data is + * present ({@code mask > 0}) vs. missing ({@code mask == 0}), and (c) the acquisition xlog zarr, which provides a + * per-slab point cloud of beam positions ({@code x_reference}, {@code y_reference} in microns) together with each + * point's distance to the region of interest ({@code distance_roi}). + *

+ * The client processes the full-resolution ({@code s0}) blocks of the tissue in parallel. For each block it first + * applies the cheap ROI-distance filter (interpolating {@code distance_roi} at the block center via inverse-distance + * weighting; blocks whose interpolated distance is {@code >= maxRoiDistance} are skipped without any I/O), then reads + * the block (skipping absent/empty blocks) and fills every pixel where the mask is 0 with the average of the tissue in + * the sections above and below (copying the single neighbor at the volume's z-boundary). Before overwriting a modified + * block, the original block is copied verbatim into a sibling {@code _backup} N5 container. Finally, only the pyramid + * blocks affected by the modified {@code s0} blocks are re-downsampled (again backing up the originals first), using + * the same per-block averaging as {@code N5DownsamplerSpark} so the pyramid stays consistent. */ public class Wafer6061Inpainter { @@ -49,354 +71,661 @@ public class Wafer6061Inpainter { public static class Parameters extends CommandLineParameters { @Parameter( names = "--n5Path", - description = "Path to the N5 container containing the data and the mask.", + description = "Path to the N5 container holding the tissue and mask (local path or gs://...).", required = true) public String n5Path; @Parameter( names = "--dataset", - description = "Name of the dataset to inpaint; assumed to be a multiscale pyramid, only s0 is inpainted.", + description = "Name of the tissue multiscale group; the full-resolution data is at /s0.", required = true) public String dataset; @Parameter( names = "--mask", - description = "Name of the mask dataset. This is supposed be a binary uint8 mask covering the whole dataset.", + description = "Name of the binary uint8 mask dataset (same grid as /s0). Pass /s0 if the mask is itself a pyramid.", required = true) public String mask; @Parameter( - names = "--output", - description = "Name of the dataset to write the inpainted data to. Only blocks that are inpainted are written. " - + "If omitted, input blocks are overwritten.") - public String output; + names = "--xlogPath", + description = "Path to the acquisition xlog zarr, e.g. xlog_wafer_61.zarr (local path or gs://...).", + required = true) + public String xlogPath; @Parameter( - names = "--inpaintingSize", - description = "Rough size of the inpainting region in pixels. This is used to determine the regions to inpaint, so better to be too large than too small.", + names = "--serial", + description = "Serial label (id_serial) of the section stored in this N5; resolved to the reference arrays' slab position.", required = true) - public int stepSize; + public long serial; + @Parameter( + names = "--scale", + description = "Nanometers per pixel used to map a block center to the xlog micron frame: micron = (pixel + translate) * scale / 1000.") + public double scale = 8.0; - public void validate() { - if (stepSize <= 0) { - throw new IllegalArgumentException("Inpainting size must be positive"); - } - } + @Parameter( + names = "--maxRoiDistance", + description = "Keep (inpaint) a block only if its interpolated distance_roi (microns) is less than this.") + public double maxRoiDistance = 10.0; + + @Parameter( + names = "--k", + description = "Number of nearest neighbors used for inverse-distance weighting.") + public int k = 8; + + @Parameter( + names = "--idwPower", + description = "Power p in the inverse-distance weights 1 / r^p.") + public double idwPower = 2.0; + + @Parameter( + names = "--downsampleFactors", + description = "Relative per-level downsampling factors of the existing pyramid, e.g. 2,2,1.") + public String downsampleFactors = "2,2,1"; + + @Parameter( + names = "--backupPath", + description = "N5 container where original (overwritten) blocks are backed up. Defaults to a sibling _backup.n5.") + public String backupPath; + + @Parameter( + names = "--dryRun", + description = "Only compute and log candidate / would-be-modified counts; do not write or back up anything.") + public boolean dryRun = false; public String fullDataset() { return dataset + "/s0"; } + + public int[] getDownsampleFactors() { + final String[] parts = downsampleFactors.split(","); + final int[] factors = new int[parts.length]; + for (int i = 0; i < parts.length; i++) { + factors[i] = Integer.parseInt(parts[i].trim()); + } + return factors; + } + + public String getBackupPath() { + if (backupPath != null) { + return backupPath; + } + String p = n5Path; + while (p.endsWith("/")) { + p = p.substring(0, p.length() - 1); + } + if (p.endsWith(".n5")) { + return p.substring(0, p.length() - 3) + "_backup.n5"; + } + return p + "_backup.n5"; + } + + public void validate() { + if (scale <= 0) { + throw new IllegalArgumentException("--scale must be positive"); + } + if (maxRoiDistance <= 0) { + throw new IllegalArgumentException("--maxRoiDistance must be positive"); + } + if (k < 1) { + throw new IllegalArgumentException("--k must be at least 1"); + } + } } - private final Parameters param; + private final Parameters params; - private ExtendedAttributes tissueAttributes; - private ExtendedAttributes maskAttributes; + public Wafer6061Inpainter(final Parameters params) { + this.params = params; + } + public static void main(final String[] args) { + final ClientRunner clientRunner = new ClientRunner(args) { + @Override + public void runClient(final String[] args) throws Exception { - public Wafer6061Inpainter(final Parameters parameters) { - this.param = parameters; + final Parameters parameters = new Parameters(); + parameters.parse(args); + parameters.validate(); + + LOG.info("runClient: entry, parameters={}", parameters); + + final Wafer6061Inpainter client = new Wafer6061Inpainter(parameters); + client.run(); + } + }; + clientRunner.run(); } - public void run() { - final String output = param.output == null ? "input dataset" : "'" + param.output + "'"; - LOG.info("Inpainting dataset '{}' in '{}' using mask '{}' and writing to {}", - param.dataset, param.n5Path, param.mask, output); - - // Read and cache some metadata of the tissue and mask datasets - // Assume that the tissue is a multiscale pyramid / mask is a standalone dataset - try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, param.n5Path)) { - LOG.info("Reading metadata from {}", param.n5Path); - tissueAttributes = ExtendedAttributes.read(n5, param.fullDataset(), param.dataset); - maskAttributes = ExtendedAttributes.read(n5, param.mask, param.mask); - - if (param.output == null) { - param.output = param.fullDataset(); - LOG.info("Output dataset equals input dataset. Overwriting blocks in the input dataset '{}'", param.output); - } else if (n5.exists(param.output)) { - throw new IllegalArgumentException("Dataset '" + param.output + "' is different from the input dataset and already exists. Stopping."); - } else { - LOG.info("Output dataset is '{}'. Creating new dataset.", param.output); - try (final N5Writer n5Writer = new N5Factory().openWriter(N5Factory.StorageFormat.N5, param.n5Path)) { - n5Writer.createDataset(param.output, tissueAttributes.attrs); + public void run() throws IOException { + + // 1. Load the (small) per-slab 2-D point cloud from the xlog on the driver. + final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial); + LOG.info("run: loaded {} ROI reference points for serial {}", cloud.size(), params.serial); + + // Read the tissue s0 metadata and discover the existing pyramid levels. + final long[] tissueTranslate; + final long[] maskTranslate; + final int numDimensions; + final List levels = new ArrayList<>(); + final Map levelAttributes = new LinkedHashMap<>(); + try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path)) { + + final DatasetAttributes s0Attributes = n5.getDatasetAttributes(params.fullDataset()); + if (s0Attributes == null) { + throw new IllegalArgumentException("tissue dataset " + params.fullDataset() + " does not exist"); + } + numDimensions = s0Attributes.getNumDimensions(); + if (numDimensions != 3) { + throw new IllegalArgumentException("expected a 3D tissue volume but " + params.fullDataset() + + " has " + numDimensions + " dimensions"); + } + if (n5.getDatasetAttributes(params.mask) == null) { + throw new IllegalArgumentException("mask dataset " + params.mask + " does not exist"); + } + + tissueTranslate = readTranslate(n5, params.fullDataset(), numDimensions); + maskTranslate = readTranslate(n5, params.mask, numDimensions); + + levels.add(params.fullDataset()); + levelAttributes.put(params.fullDataset(), s0Attributes); + for (int scale = 1; ; scale++) { + final String levelDataset = params.dataset + "/s" + scale; + if (! n5.datasetExists(levelDataset)) { + break; + } + levels.add(levelDataset); + levelAttributes.put(levelDataset, n5.getDatasetAttributes(levelDataset)); + } + } + LOG.info("run: tissue translate={}, mask translate={}, pyramid levels={}", + Arrays.toString(tissueTranslate), Arrays.toString(maskTranslate), levels); + + // Create the backup container and mirror all datasets that might receive backups. + final String backupPath = params.getBackupPath(); + if (! params.dryRun) { + try (final N5Writer backup = new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + for (final String levelDataset : levels) { + if (! backup.datasetExists(levelDataset)) { + backup.createDataset(levelDataset, levelAttributes.get(levelDataset)); + } } } + LOG.info("run: originals will be backed up to {}", backupPath); } - final SparkConf conf = new SparkConf().setAppName("Wafer6061Inpainter"); + final SparkConf conf = new SparkConf().setAppName(getClass().getSimpleName()); try (final JavaSparkContext sparkContext = new JavaSparkContext(conf)) { - runWithSparkContext(sparkContext); + LOG.info("run: appId is {}", sparkContext.getConf().getAppId()); + runWithSparkContext(sparkContext, cloud, tissueTranslate, maskTranslate, + levels, levelAttributes, backupPath); } } - private void runWithSparkContext(final JavaSparkContext sparkContext) { - // Find out which blocks need inpainting (i.e., find blocks that are neither all mask nor all void) - final List maskBlocks = Grid.create(maskAttributes.attrs.getDimensions(), maskAttributes.attrs.getBlockSize()); - final JavaRDD maskRDD = sparkContext.parallelize(maskBlocks); - final Broadcast maskAttributesBroadcast = sparkContext.broadcast(maskAttributes); - final Broadcast paramBroadcast = sparkContext.broadcast(param); - LOG.info("Filtering empty mask blocks from {} blocks", maskBlocks.size()); - - final List nonHomogeneousMaskBlocks = maskRDD - .map(block -> translateAndCheckHomogeneity(block, - maskAttributesBroadcast.value().min, - paramBroadcast.value())) - .filter(Objects::nonNull) - .collect(); - LOG.info("Found {} non-homogeneous mask blocks", nonHomogeneousMaskBlocks.size()); - - // Check which tissue blocks are covered by the potentially inpainted mask blocks determined above - final List tissueBlocks = Grid.create(tissueAttributes.attrs.getDimensions(), tissueAttributes.attrs.getBlockSize()); - final JavaRDD tissueBlocksRDD = sparkContext.parallelize(tissueBlocks); - final Broadcast tissueAttributesBroadcast = sparkContext.broadcast(tissueAttributes); - final Broadcast> maskBlocksBroadcast = sparkContext.broadcast(nonHomogeneousMaskBlocks); - - final List tissueBlocksToInpaint = tissueBlocksRDD.map( - block -> translateAndCheckIfOverlaps(block, - tissueAttributesBroadcast.value().min, - maskBlocksBroadcast.value())) - .filter(Objects::nonNull) - .collect(); - LOG.info("Found {} tissue blocks to inpaint", tissueBlocksToInpaint.size()); - - // Inpaint the blocks - final JavaRDD inpaintingBlocksRDD = sparkContext.parallelize(tissueBlocksToInpaint); - - inpaintingBlocksRDD.foreach(block -> inpaintBlock(block, - maskAttributesBroadcast.value().min, - tissueAttributesBroadcast.value().min, - paramBroadcast.value(), - tissueAttributesBroadcast.value().attrs)); + private void runWithSparkContext(final JavaSparkContext sparkContext, + final PointCloud cloud, + final long[] tissueTranslate, + final long[] maskTranslate, + final List levels, + final Map levelAttributes, + final String backupPath) { + + final DatasetAttributes s0Attributes = levelAttributes.get(params.fullDataset()); + final List s0Blocks = Grid.create(s0Attributes.getDimensions(), s0Attributes.getBlockSize()); + LOG.info("runWithSparkContext: {} s0 grid blocks to consider", s0Blocks.size()); + + final Broadcast cloudBroadcast = sparkContext.broadcast(cloud); + final Broadcast paramsBroadcast = sparkContext.broadcast(params); + + // 2. + 3. Filter (distance, then presence) and inpaint in one distributed pass. + final String backup = params.dryRun ? null : backupPath; + final JavaRDD modifiedRDD = sparkContext.parallelize(s0Blocks).mapPartitions( + blockIterator -> inpaintPartition(blockIterator, + cloudBroadcast.value(), + paramsBroadcast.value(), + tissueTranslate, + maskTranslate, + backup)); + + final List modifiedS0 = modifiedRDD.collect(); + LOG.info("runWithSparkContext: {} s0 block(s) were {}", + modifiedS0.size(), params.dryRun ? "identified for inpainting (dry run)" : "inpainted"); + + if (params.dryRun || modifiedS0.isEmpty()) { + return; + } + + // 6. Selectively update the downsample pyramid, one level at a time. + updatePyramid(sparkContext, levels, levelAttributes, modifiedS0, backupPath); } - private static Grid.Block translateAndCheckHomogeneity( - final Grid.Block block, - final long[] shift, - final Parameters param - ) { - LogUtilities.setupExecutorLog4j(""); - - // Read the mask block and check if it is homogeneous - boolean isHomogeneous = true; - try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, param.n5Path)) { - final Img mask = N5Utils.open(n5, param.mask); - final Interval interval = Intervals.intersect(mask, block); - final RandomAccessibleInterval maskPixels = Views.interval(mask, interval); - - final UnsignedByteType firstPixel = maskPixels.firstElement(); - for (final UnsignedByteType pixel : maskPixels) { - if (! pixel.equals(firstPixel)) { - isHomogeneous = false; - break; + // ------------------------------------------------------------------------------------------------ + // Step 1: load the per-slab 2-D point cloud from the xlog zarr. + // ------------------------------------------------------------------------------------------------ + + /** + * Reads {@code x_reference}, {@code y_reference} and {@code distance_roi} for the slab whose {@code id_serial} + * label equals {@code serial}, flattening the mfov x sfov grid into a single 2-D cloud (NaN points dropped). + */ + static PointCloud loadPointCloud(final String xlogPath, final long serial) { + try (final N5Reader xlog = new N5Factory().openReader(xlogPath)) { + + final double[] idSerial = read1d(xlog, "id_serial"); + final int slabPosition = findSlabPosition(idSerial, serial); + if (slabPosition < 0) { + final long[] available = new long[idSerial.length]; + for (int i = 0; i < idSerial.length; i++) { + available[i] = Math.round(idSerial[i]); } + throw new IllegalArgumentException("serial " + serial + " not found in id_serial; available serials are " + + Arrays.toString(available)); } - } - - // Translate the block to physical coordinates - final Interval blockInterval = Intervals.translate(block, shift); - final Grid.Block translatedBlock = new Grid.Block(blockInterval, block.gridPosition); - final String blockType = isHomogeneous ? "homogeneous -> skip" : "non-homogeneous -> possibly inpaint"; - LOG.info("Mask block {} at {} is {}", translatedBlock.gridPosition, translatedBlock.offset, blockType); - return isHomogeneous ? null : translatedBlock; - } + // The slab axis is the one whose size matches the length of id_serial (413 for wafer 61). + final long slabCount = idSerial.length; + final RandomAccessibleInterval xSlab = readSlab(xlog, "x_reference", slabCount, slabPosition); + final RandomAccessibleInterval ySlab = readSlab(xlog, "y_reference", slabCount, slabPosition); + final RandomAccessibleInterval distSlab = readSlab(xlog, "distance_roi", slabCount, slabPosition); + + final List xs = new ArrayList<>(); + final List ys = new ArrayList<>(); + final List dists = new ArrayList<>(); + + final Cursor xc = xSlab.localizingCursor(); + final RandomAccess yra = ySlab.randomAccess(); + final RandomAccess dra = distSlab.randomAccess(); + final long[] pos = new long[xSlab.numDimensions()]; + while (xc.hasNext()) { + final double x = xc.next().get(); + xc.localize(pos); + final double y = yra.setPositionAndGet(pos).get(); + final double d = dra.setPositionAndGet(pos).get(); + if (Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(d)) { + xs.add(x); + ys.add(y); + dists.add(d); + } + } - private static Grid.Block translateAndCheckIfOverlaps( - final Grid.Block block, - final long[] shift, - final List blocksToCheckAgainst - ) { - LogUtilities.setupExecutorLog4j(""); - - // Translate the block to physical coordinates - final Interval blockInterval = Intervals.translate(block, shift); - final Grid.Block translatedBlock = new Grid.Block(blockInterval, block.gridPosition); - - // Check if the block overlaps with any of the mask blocks that might need inpainting - for (final Interval maskBlock : blocksToCheckAgainst) { - final boolean intervalsAreDisjoint = Intervals.isEmpty(Intervals.intersect(translatedBlock, maskBlock)); - if (! intervalsAreDisjoint) { - LOG.info("Tissue block {} at {} is determined a candidate for inpainting", - translatedBlock.gridPosition, translatedBlock.minAsLongArray()); - return translatedBlock; + if (xs.isEmpty()) { + throw new IllegalArgumentException("no finite ROI reference points found for serial " + serial + + " (slab position " + slabPosition + ")"); } + return new PointCloud(toArray(xs), toArray(ys), toArray(dists)); } + } - LOG.info("Tissue block {} at {} is not a candidate for inpainting", - translatedBlock.gridPosition, translatedBlock.minAsLongArray()); - return null; + /** Reads a 1-D double dataset (blosc-safe: uses the block-not-found overload to avoid the getAttribute NPE). */ + private static double[] read1d(final N5Reader n5, final String dataset) { + final RandomAccessibleInterval img = openDoubles(n5, dataset); + final long length = img.dimension(0); + final double[] values = new double[(int) length]; + final RandomAccess ra = img.randomAccess(); + for (int i = 0; i < length; i++) { + values[i] = ra.setPositionAndGet(new long[] {i}).get(); + } + return values; } - private static void inpaintBlock( - final Grid.Block block, - final long[] maskMin, - final long[] tissueMin, - final Parameters param, - final DatasetAttributes targetAttributes - ) { - LogUtilities.setupExecutorLog4j("Block " + Arrays.toString(block.gridPosition)); - - // Preallocate the inpainted block - final Img inpaintedBlock = ArrayImgs.unsignedBytes(block.dimensions); - - try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, param.n5Path)) { - // Load and translate the tissue and mask data - LOG.info("Loading data at {}", block.offset); - final Img rawTissue = N5Utils.open(n5, param.fullDataset()); - final Img rawMask = N5Utils.open(n5, param.mask); - - final RandomAccessibleInterval tissue = Views.translate(rawTissue, tissueMin); - final RandomAccessibleInterval mask = Views.translate(rawMask, maskMin); - - // For each pixel, determine if it should be inpainted and if so, inpaint it by interpolating in z - LOG.info("Start inpainting"); - final long start = System.currentTimeMillis(); - final Cursor targetCursor = Views.translate(inpaintedBlock, block.offset).localizingCursor(); - final long[] location = new long[3]; - final PixelFiller interpolator = new PixelFiller(tissue, mask, param.stepSize); - - while (targetCursor.hasNext()) { - final UnsignedByteType targetPixel = targetCursor.next(); - targetCursor.localize(location); - final int value = interpolator.getPixel(location); - targetPixel.set(value); + /** Opens a 3-D reference array and returns the 2-D slice for the given slab position. */ + private static RandomAccessibleInterval readSlab(final N5Reader n5, + final String dataset, + final long slabCount, + final int slabPosition) { + final RandomAccessibleInterval img = openDoubles(n5, dataset); + int slabAxis = -1; + for (int d = 0; d < img.numDimensions(); d++) { + if (img.dimension(d) == slabCount) { + slabAxis = d; + break; } - LOG.info("Finished inpainting in {} ms", System.currentTimeMillis() - start); } - - try (final N5Writer n5Writer = new N5Factory().openWriter(N5Factory.StorageFormat.N5, param.n5Path)) { - N5Utils.saveBlock(inpaintedBlock, n5Writer, param.output, targetAttributes, block.gridPosition); - LOG.info("Wrote tissue block to '{}'", param.output); - } catch (final Exception e) { - LOG.error("Failed to write inpainted block", e); + if (slabAxis < 0) { + throw new IllegalArgumentException("could not find slab axis (size " + slabCount + ") in " + dataset + + " with dimensions " + Arrays.toString(img.dimensionsAsLongArray())); } + return Views.hyperSlice(img, slabAxis, slabPosition); } + /** + * Opens a double-typed dataset using the {@code (blockNotFoundHandler, accessFlags)} overload. This avoids + * {@code N5Utils.open(reader, dataset)}'s {@code isLabelMultisetType -> getAttribute} path, which NPEs on blosc + * arrays in the shaded jar (the n5-blosc CompressionType service registration is filtered out). Missing blocks + * are filled with NaN so they are dropped downstream. + */ + private static RandomAccessibleInterval openDoubles(final N5Reader n5, final String dataset) { + final Consumer> nanFill = it -> it.forEach(t -> t.set(Double.NaN)); + return N5Utils.open(n5, dataset, nanFill, AccessFlags.setOf()); + } - public static void main(final String[] args) { - final ClientRunner clientRunner = new ClientRunner(args) { - @Override - public void runClient(final String[] args) { + // ------------------------------------------------------------------------------------------------ + // Steps 2 + 3: filter and inpaint one partition of s0 blocks. + // ------------------------------------------------------------------------------------------------ - final Wafer6061Inpainter.Parameters parameters = new Wafer6061Inpainter.Parameters(); - parameters.parse(args); - parameters.validate(); + private static Iterator inpaintPartition(final Iterator blocks, + final PointCloud cloud, + final Parameters params, + final long[] tissueTranslate, + final long[] maskTranslate, + final String backupPath) { - LOG.info("runClient: entry, parameters={}", parameters); + LogUtilities.setupExecutorLog4j("inpaint"); + + final List modified = new ArrayList<>(); + final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(params.k); + final double micronPerPixel = params.scale / 1000.0; + + try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); + final N5Writer tissueWriter = params.dryRun ? null : + new N5Factory().openWriter(N5Factory.StorageFormat.N5, params.n5Path); + final N5Writer backupWriter = params.dryRun ? null : + new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + + final DatasetAttributes s0Attributes = reader.getDatasetAttributes(params.fullDataset()); + + final Img rawTissue = N5Utils.open(reader, params.fullDataset()); + final Img rawMask = N5Utils.open(reader, params.mask); + final RandomAccessibleInterval tissue = Views.translate(rawTissue, tissueTranslate); + final RandomAccessibleInterval mask = Views.translate(rawMask, maskTranslate); - final Wafer6061Inpainter inpainter = new Wafer6061Inpainter(parameters); - inpainter.run(); + final long zMin = tissueTranslate[2]; + final long zMax = tissueTranslate[2] + s0Attributes.getDimensions()[2] - 1; + + while (blocks.hasNext()) { + final Grid.Block block = blocks.next(); + + // (1) ROI-distance filter, computed with no I/O. + final double centerX = (block.offset[0] + block.dimensions[0] / 2.0 + tissueTranslate[0]) * micronPerPixel; + final double centerY = (block.offset[1] + block.dimensions[1] / 2.0 + tissueTranslate[1]) * micronPerPixel; + final double roiDistance = cloud.interpolate(search, centerX, centerY, params.idwPower); + if (! (roiDistance < params.maxRoiDistance)) { + continue; + } + + // (2) presence check: readBlock returns null for absent (empty) blocks. + final DataBlock originalBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); + if (originalBlock == null) { + continue; + } + + // (3) inpaint: fill mask==0 pixels with the z-average of the sections above/below. + final InpaintResult result = inpaintBlock(tissue, mask, block, tissueTranslate, zMin, zMax); + if (! result.changed) { + continue; + } + + modified.add(block.gridPosition); + if (! params.dryRun) { + backupWriter.writeBlock(params.fullDataset(), s0Attributes, originalBlock); + N5Utils.saveBlock(result.inpainted, tissueWriter, params.fullDataset(), s0Attributes, block.gridPosition); + } + LOG.info("inpaintPartition: inpainted block {} (roiDistance={})", + Arrays.toString(block.gridPosition), roiDistance); } - }; - clientRunner.run(); - } + } + return modified.iterator(); + } /** - * Performs all the inpainting-logic, i.e., when and how to inpaint. + * Produces the inpainted version of a single block: every pixel where the mask is 0 is replaced with the average + * of the tissue in the sections above and below (the single neighbor is copied at the volume's z-boundary). The + * {@code tissue} and {@code mask} are expected to be world-translated (i.e. indexed in the same coordinate frame, + * offset by their respective {@code translate}); {@code zMin}/{@code zMax} are the world z-bounds of the tissue. */ - private static class PixelFiller { - - private final RandomAccess tissueAccess; - private final RandomAccess maskAccess; - private final int posStep; - private final int negStep; - - public PixelFiller( - final RandomAccessibleInterval tissue, - final RandomAccessibleInterval mask, - final int stepSize - ) { - this.tissueAccess = tissue.randomAccess(); - this.maskAccess = Views.extendZero(mask).randomAccess(); - this.posStep = stepSize; - this.negStep = -2 * stepSize; - } - - public int getPixel(final long[] position) { - if (shouldBeInpainted(position)) { - return zAverage(position); + static InpaintResult inpaintBlock(final RandomAccessibleInterval tissue, + final RandomAccessibleInterval mask, + final Grid.Block block, + final long[] tissueTranslate, + final long zMin, + final long zMax) { + + final RandomAccess tissueAccess = tissue.randomAccess(); + final RandomAccess maskAccess = Views.extendZero(mask).randomAccess(); + + final Img inpainted = ArrayImgs.unsignedBytes(block.dimensions); + final Cursor cursor = inpainted.localizingCursor(); + final long[] local = new long[3]; + final long[] world = new long[3]; + boolean changed = false; + while (cursor.hasNext()) { + final UnsignedByteType target = cursor.next(); + cursor.localize(local); + world[0] = block.offset[0] + tissueTranslate[0] + local[0]; + world[1] = block.offset[1] + tissueTranslate[1] + local[1]; + world[2] = block.offset[2] + tissueTranslate[2] + local[2]; + + final int original = tissueAccess.setPositionAndGet(world).get(); + final int value; + if (maskAccess.setPositionAndGet(world).get() > 0) { + value = original; } else { - return tissueAccess.setPositionAndGet(position).get(); + value = zAverage(tissueAccess, world, zMin, zMax); + if (value != original) { + changed = true; + } } + target.set(value); } - /** - * Average the z-values of the pixels above and below the current pixel. - * If only one of the pixels is available, that value is used. - * If neither is available, the pixel is set to 0. - */ - private int zAverage(final long[] position) { - maskAccess.setPosition(position); - maskAccess.move(-1, 2); - final boolean hasContentAbove = maskAccess.get().get() > 0; - maskAccess.move(2, 2); - final boolean hasContentBelow = maskAccess.get().get() > 0; - - tissueAccess.setPositionAndGet(position); - if (hasContentAbove && hasContentBelow) { - tissueAccess.move(-1, 2); - final int above = tissueAccess.get().get(); - - tissueAccess.move(2, 2); - final int below = tissueAccess.get().get(); - - return UnsignedByteType.getCodedSignedByteChecked((above + below) >>> 1); - } else if (hasContentAbove) { - tissueAccess.move(-1, 2); - return tissueAccess.get().get(); - } else if (hasContentBelow) { - tissueAccess.move(2, 2); - return tissueAccess.get().get(); - } else { - return 0; + return new InpaintResult(inpainted, changed); + } + + /** Average of the tissue in the z-1 and z+1 sections, copying the single neighbor at the volume z-boundary. */ + private static int zAverage(final RandomAccess tissueAccess, + final long[] world, + final long zMin, + final long zMax) { + final long z = world[2]; + final boolean hasAbove = (z - 1) >= zMin; + final boolean hasBelow = (z + 1) <= zMax; + if (hasAbove && hasBelow) { + world[2] = z - 1; + final int above = tissueAccess.setPositionAndGet(world).get(); + world[2] = z + 1; + final int below = tissueAccess.setPositionAndGet(world).get(); + world[2] = z; + return (above + below) >>> 1; + } else if (hasAbove) { + world[2] = z - 1; + final int above = tissueAccess.setPositionAndGet(world).get(); + world[2] = z; + return above; + } else if (hasBelow) { + world[2] = z + 1; + final int below = tissueAccess.setPositionAndGet(world).get(); + world[2] = z; + return below; + } else { + return tissueAccess.setPositionAndGet(world).get(); + } + } + + // ------------------------------------------------------------------------------------------------ + // Step 6: selectively re-downsample only the pyramid blocks affected by the modified s0 blocks. + // ------------------------------------------------------------------------------------------------ + + private void updatePyramid(final JavaSparkContext sparkContext, + final List levels, + final Map levelAttributes, + final List modifiedS0, + final String backupPath) { + + final int[] factors = params.getDownsampleFactors(); + List modifiedPrevious = modifiedS0; + + for (int scale = 1; scale < levels.size(); scale++) { + final String fromDataset = levels.get(scale - 1); + final String toDataset = levels.get(scale); + final DatasetAttributes toAttributes = levelAttributes.get(toDataset); + + // affected blocks: previous-level grid position p maps to this-level block p / factor. + final Map affected = new LinkedHashMap<>(); + for (final long[] p : modifiedPrevious) { + final long[] g = affectedBlock(p, factors); + affected.putIfAbsent(Arrays.toString(g), g); } + final List affectedBlocks = new ArrayList<>(affected.values()); + LOG.info("updatePyramid: re-downsampling {} block(s) for {}", affectedBlocks.size(), toDataset); + + final String n5Path = params.n5Path; + sparkContext.parallelize(affectedBlocks).foreach( + gridPosition -> downsampleBlock(gridPosition, n5Path, backupPath, + fromDataset, toDataset, toAttributes, factors)); + + modifiedPrevious = affectedBlocks; + } + } + + /** + * Re-downsamples a single pyramid block from the (already updated) previous level, replicating the per-block + * math of {@code N5DownsamplerSpark} so the result matches the rest of the pyramid. The original block is backed + * up before being overwritten. + */ + static void downsampleBlock(final long[] gridPosition, + final String n5Path, + final String backupPath, + final String fromDataset, + final String toDataset, + final DatasetAttributes toAttributes, + final int[] factors) { + + LogUtilities.setupExecutorLog4j("downsample"); + + final int n = toAttributes.getNumDimensions(); + final CellGrid cellGrid = new CellGrid(toAttributes.getDimensions(), toAttributes.getBlockSize()); + + final long[] targetMin = new long[n]; + final int[] cellDimensions = new int[n]; + cellGrid.getCellDimensions(gridPosition, targetMin, cellDimensions); + + final long[] sourceMin = new long[n]; + final long[] sourceSize = new long[n]; + final long[] targetSize = new long[n]; + for (int d = 0; d < n; d++) { + sourceMin[d] = targetMin[d] * factors[d]; + sourceSize[d] = (long) cellDimensions[d] * factors[d]; + targetSize[d] = cellDimensions[d]; } - /** - * Determines if the pixel at the given position should be inpainted based on the local environment. - */ - private boolean shouldBeInpainted(final long[] position) { - final boolean hasContent = maskAccess.setPositionAndGet(position).get() > 0; - if (hasContent) { - return false; + try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, n5Path); + final N5Writer writer = new N5Factory().openWriter(N5Factory.StorageFormat.N5, n5Path); + final N5Writer backupWriter = new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + + final RandomAccessibleInterval source = N5Utils.open(reader, fromDataset); + final RandomAccessibleInterval sourceBlock = Views.offsetInterval(source, sourceMin, sourceSize); + + final Img targetBlock = ArrayImgs.unsignedBytes(targetSize); + Downsample.downsample(sourceBlock, targetBlock, factors); + + // back up the original block (if present) before overwriting. + final DataBlock originalBlock = reader.readBlock(toDataset, toAttributes, gridPosition); + if (originalBlock != null) { + backupWriter.writeBlock(toDataset, toAttributes, originalBlock); } - // If the pixel has no content, check the pixels in +/- y direction - // Only if both have content, the pixel should be inpainted (otherwise, it is a border pixel) - maskAccess.move(posStep, 1); - final boolean hasContentFront = maskAccess.get().get() > 0; - maskAccess.move(negStep, 1); - final boolean hasContentBack = maskAccess.get().get() > 0; - if (hasContentFront && hasContentBack) { - return true; + // delete first so a block that became empty does not leave a stale remnant. + N5Utils.deleteBlock(targetBlock, writer, toDataset, gridPosition); + N5Utils.saveNonEmptyBlock(targetBlock, writer, toDataset, gridPosition, new UnsignedByteType()); + LOG.info("downsampleBlock: updated {} block {}", toDataset, Arrays.toString(gridPosition)); + } + } + + // ------------------------------------------------------------------------------------------------ + // Helpers + // ------------------------------------------------------------------------------------------------ + + /** Returns the slab-axis position whose id_serial label equals {@code serial}, or -1 if none matches. */ + static int findSlabPosition(final double[] idSerial, final long serial) { + for (int i = 0; i < idSerial.length; i++) { + if (Math.round(idSerial[i]) == serial) { + return i; } + } + return -1; + } - // If that is inconclusive, check the pixels in +/- x direction - maskAccess.move(posStep, 1); - maskAccess.move(posStep, 0); - final boolean hasContentRight = maskAccess.get().get() > 0; - maskAccess.move(negStep, 0); - final boolean hasContentLeft = maskAccess.get().get() > 0; - return hasContentRight && hasContentLeft; + /** Maps a grid block position at pyramid level k-1 to the block it feeds at level k (per-dimension p / factor). */ + static long[] affectedBlock(final long[] previousGridPosition, final int[] factors) { + final long[] g = new long[previousGridPosition.length]; + for (int d = 0; d < g.length; d++) { + g[d] = previousGridPosition[d] / factors[d]; } + return g; + } + + private static long[] readTranslate(final N5Reader n5, final String dataset, final int numDimensions) { + final long[] translate = n5.getAttribute(dataset, "translate", long[].class); + return translate != null ? translate : new long[numDimensions]; + } + + private static double[] toArray(final List values) { + final double[] array = new double[values.size()]; + for (int i = 0; i < array.length; i++) { + array[i] = values.get(i); + } + return array; + } + /** Result of inpainting one block: the (block-sized) inpainted image and whether any pixel changed. */ + static class InpaintResult { + final Img inpainted; + final boolean changed; + + InpaintResult(final Img inpainted, final boolean changed) { + this.inpainted = inpainted; + this.changed = changed; + } } - private static class ExtendedAttributes implements Serializable { - public final DatasetAttributes attrs; - public final long[] min; + /** Serializable 2-D point cloud of ROI reference points with a distance value per point. */ + static class PointCloud implements Serializable { + private final double[] xs; + private final double[] ys; + private final double[] dists; - public ExtendedAttributes(final DatasetAttributes attrs, final long[] min) { - this.attrs = attrs; - this.min = min; + PointCloud(final double[] xs, final double[] ys, final double[] dists) { + this.xs = xs; + this.ys = ys; + this.dists = dists; } - public static ExtendedAttributes read(final N5Reader n5, final String attrsPath, final String minPath) { - final DatasetAttributes attrs = n5.getDatasetAttributes(attrsPath); - final long[] min = n5.getAttribute(minPath, "translate", long[].class); - return new ExtendedAttributes(attrs, min); + int size() { + return xs.length; + } + + KNearestNeighborSearchOnKDTree buildSearch(final int k) { + final List points = new ArrayList<>(xs.length); + final List values = new ArrayList<>(xs.length); + for (int i = 0; i < xs.length; i++) { + points.add(new RealPoint(xs[i], ys[i])); + values.add(new DoubleType(dists[i])); + } + final KDTree tree = new KDTree<>(values, points); + return new KNearestNeighborSearchOnKDTree<>(tree, Math.min(k, xs.length)); + } + + /** Inverse-distance-weighted interpolation of the distance value at (x, y). */ + double interpolate(final KNearestNeighborSearchOnKDTree search, + final double x, + final double y, + final double power) { + search.search(new RealPoint(x, y)); + final int numNeighbors = search.getK(); + double numerator = 0; + double denominator = 0; + for (int i = 0; i < numNeighbors; i++) { + final double r = search.getDistance(i); + final double v = search.getSampler(i).get().get(); + if (r == 0.0) { + return v; + } + final double w = 1.0 / Math.pow(r, power); + numerator += w * v; + denominator += w; + } + return numerator / denominator; } } } From d355386a876e4e7ebbae5fd80685a9e7947d701a Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sat, 11 Jul 2026 16:28:08 -0400 Subject: [PATCH 02/15] Fix coordinate transformations of ROI --- .../spark/multisem/Wafer6061Inpainter.java | 241 ++++++++++++++---- 1 file changed, 194 insertions(+), 47 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index f490d66cf..4452f30e5 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -51,9 +51,11 @@ * Spark client for filling holes in the tissue of wafer 60/61 N5 volumes. *

* Inputs are (a) an N5 tissue volume (only non-empty blocks stored), (b) an N5 mask marking where image data is - * present ({@code mask > 0}) vs. missing ({@code mask == 0}), and (c) the acquisition xlog zarr, which provides a - * per-slab point cloud of beam positions ({@code x_reference}, {@code y_reference} in microns) together with each - * point's distance to the region of interest ({@code distance_roi}). + * present ({@code mask > 0}) vs. missing ({@code mask == 0}), and (c) the acquisition xlog zarr. The ROI point cloud + * is built in the tissue's s0 voxel frame directly from the xlog: each SFOV's acquisition position ({@code x}, + * {@code y}) is placed exactly as {@code msem_to_render.py} ingests it (rotate by {@code 180 + rotation_slab}, then + * {@code stage - min + halfSFOV}), and carries that SFOV's {@code distance_roi}. No render service and no coordinate + * fitting are needed (x/y/distance_roi share the same {@code (mfov, sfov)} indexing). *

* The client processes the full-resolution ({@code s0}) blocks of the tissue in parallel. For each block it first * applies the cheap ROI-distance filter (interpolating {@code distance_roi} at the block center via inverse-distance @@ -100,9 +102,20 @@ public static class Parameters extends CommandLineParameters { public long serial; @Parameter( - names = "--scale", - description = "Nanometers per pixel used to map a block center to the xlog micron frame: micron = (pixel + translate) * scale / 1000.") - public double scale = 8.0; + names = "--sfovWidth", + description = "SFOV image width in pixels (xlog X_SFOV size); used for the ingestion half-SFOV placement offset.") + public int sfovWidth = 2000; + + @Parameter( + names = "--sfovHeight", + description = "SFOV image height in pixels (xlog Y_SFOV size); used for the ingestion half-SFOV placement offset.") + public int sfovHeight = 1748; + + @Parameter( + names = "--scan", + description = "xlog scan index whose x/y positions define the cloud. Default (-1) auto-picks the first scan " + + "with finite x and rotation_slab for the slab (positions are nearly scan-independent).") + public int scan = -1; @Parameter( names = "--maxRoiDistance", @@ -162,8 +175,8 @@ public String getBackupPath() { } public void validate() { - if (scale <= 0) { - throw new IllegalArgumentException("--scale must be positive"); + if (sfovWidth <= 0 || sfovHeight <= 0) { + throw new IllegalArgumentException("--sfovWidth and --sfovHeight must be positive"); } if (maxRoiDistance <= 0) { throw new IllegalArgumentException("--maxRoiDistance must be positive"); @@ -200,8 +213,9 @@ public void runClient(final String[] args) throws Exception { public void run() throws IOException { - // 1. Load the (small) per-slab 2-D point cloud from the xlog on the driver. - final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial); + // 1. Load the (small) per-slab 2-D point cloud from the xlog on the driver, placed in the tissue voxel frame. + final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial, + params.scan, params.sfovWidth, params.sfovHeight); LOG.info("run: loaded {} ROI reference points for serial {}", cloud.size(), params.serial); // Read the tissue s0 metadata and discover the existing pyramid levels. @@ -305,10 +319,24 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, // ------------------------------------------------------------------------------------------------ /** - * Reads {@code x_reference}, {@code y_reference} and {@code distance_roi} for the slab whose {@code id_serial} - * label equals {@code serial}, flattening the mfov x sfov grid into a single 2-D cloud (NaN points dropped). + * Builds the per-slab ROI point cloud in the tissue s0 voxel frame, entirely from the xlog. For the slab + * whose {@code id_serial} equals {@code serial}, each SFOV's acquisition position ({@code x}, {@code y}, in + * full-resolution pixels) is placed into the render/ingestion frame exactly as {@code msem_to_render.py} does: + * 

+	 *   center = EuclideanTransform(rotation = radians(180 + rotation_slab)) . (x, y)
+	 *   voxel  = center - min(center) + (sfovWidth/2, sfovHeight/2)
+	 *
+ * i.e. the ingestion {@code stage - min + margin} placement (the constant {@code margin} and the export + * {@code translate} cancel out into the voxel frame). Each voxel carries the scan-independent {@code distance_roi} + * for that SFOV. Since {@code x}, {@code y} and {@code distance_roi} are all indexed by {@code (mfov, sfov)} in the + * xlog, the correspondence is exact and no fitting is needed. Alignment (montage stitching) only perturbs these + * positions slightly, well within {@code maxRoiDistance}, so the unaligned ingestion placement is used directly. */ - static PointCloud loadPointCloud(final String xlogPath, final long serial) { + static PointCloud loadPointCloud(final String xlogPath, + final long serial, + final int scanOverride, + final int sfovWidth, + final int sfovHeight) { try (final N5Reader xlog = new N5Factory().openReader(xlogPath)) { final double[] idSerial = read1d(xlog, "id_serial"); @@ -321,41 +349,159 @@ static PointCloud loadPointCloud(final String xlogPath, final long serial) { throw new IllegalArgumentException("serial " + serial + " not found in id_serial; available serials are " + Arrays.toString(available)); } + final long slabCount = idSerial.length; // 413 for wafer 61 + + // SFOV image size defines the half-SFOV placement offset; prefer the xlog (x_sfov / y_sfov), fall back to args. + final int sfW = sfovSize(xlog, "x_sfov", sfovWidth); + final int sfH = sfovSize(xlog, "y_sfov", sfovHeight); - // The slab axis is the one whose size matches the length of id_serial (413 for wafer 61). - final long slabCount = idSerial.length; - final RandomAccessibleInterval xSlab = readSlab(xlog, "x_reference", slabCount, slabPosition); - final RandomAccessibleInterval ySlab = readSlab(xlog, "y_reference", slabCount, slabPosition); + // distance_roi is scan-independent: [slab, mfov, sfov] -> 2-D (sfov, mfov) after slicing the slab axis. final RandomAccessibleInterval distSlab = readSlab(xlog, "distance_roi", slabCount, slabPosition); - final List xs = new ArrayList<>(); - final List ys = new ArrayList<>(); - final List dists = new ArrayList<>(); + // x / y are [scan, slab, mfov, sfov]; rotation_slab is [scan, slab]. Identify axes by their (unique) sizes. + final RandomAccessibleInterval xAll = openDoubles(xlog, "x"); + final RandomAccessibleInterval yAll = openDoubles(xlog, "y"); + final RandomAccessibleInterval rotAll = openDoubles(xlog, "rotation_slab"); + final int rotSlabAxis = axisOfSize(rotAll, slabCount, "rotation_slab"); + final int rotScanAxis = 1 - rotSlabAxis; + final long nScans = rotAll.dimension(rotScanAxis); + final int slabAxisX = axisOfSize(xAll, slabCount, "x"); + final int scanAxisX = axisOfSize(xAll, nScans, "x"); + + // Choose the scan whose x/y define the cloud (positions are nearly scan-independent). + int scan = scanOverride; + if (scan < 0) { + for (int s = 0; s < nScans; s++) { + if (! Double.isFinite(rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, s))) { + continue; + } + if (hasFinite(sliceScanAndSlab(xAll, scanAxisX, s, slabAxisX, slabPosition))) { + scan = s; + break; + } + } + if (scan < 0) { + throw new IllegalArgumentException("no scan with finite x and rotation_slab found for serial " + + serial + " (slab position " + slabPosition + ")"); + } + } + + final double rotationSlab = rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, scan); + final double theta = Math.toRadians(180.0 + rotationSlab); + final double cos = Math.cos(theta); + final double sin = Math.sin(theta); + + final RandomAccessibleInterval xSlab = sliceScanAndSlab(xAll, scanAxisX, scan, slabAxisX, slabPosition); + final RandomAccessibleInterval ySlab = sliceScanAndSlab(yAll, scanAxisX, scan, slabAxisX, slabPosition); - final Cursor xc = xSlab.localizingCursor(); + // Place each SFOV center (rotation only) and keep the finite ones with their distance_roi. + final List cxs = new ArrayList<>(); + final List cys = new ArrayList<>(); + final List dists = new ArrayList<>(); + final RandomAccess xra = xSlab.randomAccess(); final RandomAccess yra = ySlab.randomAccess(); final RandomAccess dra = distSlab.randomAccess(); - final long[] pos = new long[xSlab.numDimensions()]; - while (xc.hasNext()) { - final double x = xc.next().get(); - xc.localize(pos); - final double y = yra.setPositionAndGet(pos).get(); - final double d = dra.setPositionAndGet(pos).get(); - if (Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(d)) { - xs.add(x); - ys.add(y); - dists.add(d); + final long[] pos = new long[2]; + double minX = Double.POSITIVE_INFINITY; + double minY = Double.POSITIVE_INFINITY; + for (long i0 = 0; i0 < xSlab.dimension(0); i0++) { + for (long i1 = 0; i1 < xSlab.dimension(1); i1++) { + pos[0] = i0; + pos[1] = i1; + final double x = xra.setPositionAndGet(pos).get(); + final double y = yra.setPositionAndGet(pos).get(); + final double d = dra.setPositionAndGet(pos).get(); + if (Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(d)) { + final double cx = cos * x - sin * y; + final double cy = sin * x + cos * y; + cxs.add(cx); + cys.add(cy); + dists.add(d); + minX = Math.min(minX, cx); + minY = Math.min(minY, cy); + } } } + if (cxs.isEmpty()) { + throw new IllegalArgumentException("no finite ROI points found for serial " + serial + + " (slab position " + slabPosition + ", scan " + scan + ")"); + } - if (xs.isEmpty()) { - throw new IllegalArgumentException("no finite ROI reference points found for serial " + serial + - " (slab position " + slabPosition + ")"); + // Shift into the voxel frame: min(center) -> half-SFOV (so the min tile's top-left is at the origin). + final double[] xs = new double[cxs.size()]; + final double[] ys = new double[cys.size()]; + final double[] ds = new double[dists.size()]; + for (int i = 0; i < xs.length; i++) { + xs[i] = cxs.get(i) - minX + sfW / 2.0; + ys[i] = cys.get(i) - minY + sfH / 2.0; + ds[i] = dists.get(i); } - return new PointCloud(toArray(xs), toArray(ys), toArray(dists)); + LOG.info("loadPointCloud: serial {} -> slab position {}, scan {}, rotation_slab {} deg, sfov {}x{}, {} points", + serial, slabPosition, scan, rotationSlab, sfW, sfH, xs.length); + return new PointCloud(xs, ys, ds); } } + /** SFOV image size from the xlog {@code x_sfov}/{@code y_sfov} length, or {@code fallback} if that dataset is absent. */ + private static int sfovSize(final N5Reader xlog, final String dataset, final int fallback) { + try { + final DatasetAttributes attributes = xlog.getDatasetAttributes(dataset); + if (attributes != null && attributes.getNumDimensions() >= 1) { + return (int) attributes.getDimensions()[0]; + } + } catch (final Exception e) { + LOG.warn("sfovSize: could not read {} ({}), using fallback {}", dataset, e.getMessage(), fallback); + } + return fallback; + } + + /** Index of the (unique-sized) axis matching {@code size}. */ + private static int axisOfSize(final RandomAccessibleInterval img, final long size, final String name) { + for (int d = 0; d < img.numDimensions(); d++) { + if (img.dimension(d) == size) { + return d; + } + } + throw new IllegalArgumentException("no axis of size " + size + " in " + name + " with dimensions " + + Arrays.toString(img.dimensionsAsLongArray())); + } + + /** Value of the 2-D {@code rotation_slab} at the given slab position and scan. */ + private static double rotationAt(final RandomAccessibleInterval rotAll, + final int slabAxis, + final int scanAxis, + final int slabPosition, + final int scan) { + final long[] p = new long[2]; + p[slabAxis] = slabPosition; + p[scanAxis] = scan; + return rotAll.randomAccess().setPositionAndGet(p).get(); + } + + /** Slices the 4-D {@code x}/{@code y} array to the 2-D (sfov, mfov) plane for one scan and slab. */ + private static RandomAccessibleInterval sliceScanAndSlab(final RandomAccessibleInterval img, + final int scanAxis, + final long scan, + final int slabAxis, + final long slabPosition) { + // hyper-slice the higher-indexed axis first so the lower index stays valid afterwards + final int hi = Math.max(scanAxis, slabAxis); + final int lo = Math.min(scanAxis, slabAxis); + final long hiPos = (hi == scanAxis) ? scan : slabPosition; + final long loPos = (lo == scanAxis) ? scan : slabPosition; + return Views.hyperSlice(Views.hyperSlice(img, hi, hiPos), lo, loPos); + } + + /** True if the interval has at least one finite value. */ + private static boolean hasFinite(final RandomAccessibleInterval img) { + for (final DoubleType t : Views.iterable(img)) { + if (Double.isFinite(t.get())) { + return true; + } + } + return false; + } + /** Reads a 1-D double dataset (blosc-safe: uses the block-not-found overload to avoid the getAttribute NPE). */ private static double[] read1d(final N5Reader n5, final String dataset) { final RandomAccessibleInterval img = openDoubles(n5, dataset); @@ -414,7 +560,6 @@ private static Iterator inpaintPartition(final Iterator bloc final List modified = new ArrayList<>(); final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(params.k); - final double micronPerPixel = params.scale / 1000.0; try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); final N5Writer tissueWriter = params.dryRun ? null : @@ -432,22 +577,30 @@ private static Iterator inpaintPartition(final Iterator bloc final long zMin = tissueTranslate[2]; final long zMax = tissueTranslate[2] + s0Attributes.getDimensions()[2] - 1; + long considered = 0; + long nearRoi = 0; + long present = 0; while (blocks.hasNext()) { final Grid.Block block = blocks.next(); + considered++; - // (1) ROI-distance filter, computed with no I/O. - final double centerX = (block.offset[0] + block.dimensions[0] / 2.0 + tissueTranslate[0]) * micronPerPixel; - final double centerY = (block.offset[1] + block.dimensions[1] / 2.0 + tissueTranslate[1]) * micronPerPixel; + // (1) ROI-distance filter, computed with no I/O. The cloud is already in the s0 voxel frame, so the + // block center (voxel index) queries it directly; distance_roi is interpolated by inverse-distance + // weighting (the interpolated value is in microns regardless of the voxel-space query units). + final double centerX = block.offset[0] + block.dimensions[0] / 2.0; + final double centerY = block.offset[1] + block.dimensions[1] / 2.0; final double roiDistance = cloud.interpolate(search, centerX, centerY, params.idwPower); if (! (roiDistance < params.maxRoiDistance)) { continue; } + nearRoi++; // (2) presence check: readBlock returns null for absent (empty) blocks. final DataBlock originalBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); if (originalBlock == null) { continue; } + present++; // (3) inpaint: fill mask==0 pixels with the z-average of the sections above/below. final InpaintResult result = inpaintBlock(tissue, mask, block, tissueTranslate, zMin, zMax); @@ -463,6 +616,8 @@ private static Iterator inpaintPartition(final Iterator bloc LOG.info("inpaintPartition: inpainted block {} (roiDistance={})", Arrays.toString(block.gridPosition), roiDistance); } + LOG.info("inpaintPartition: partition summary: considered={}, nearRoi(<{}um)={}, present={}, withHoles(modified)={}", + considered, params.maxRoiDistance, nearRoi, present, modified.size()); } return modified.iterator(); @@ -660,14 +815,6 @@ private static long[] readTranslate(final N5Reader n5, final String dataset, fin return translate != null ? translate : new long[numDimensions]; } - private static double[] toArray(final List values) { - final double[] array = new double[values.size()]; - for (int i = 0; i < array.length; i++) { - array[i] = values.get(i); - } - return array; - } - /** Result of inpainting one block: the (block-sized) inpainted image and whether any pixel changed. */ static class InpaintResult { final Img inpainted; From c7601a43868e28b353f07f2a73ccc651048ea091 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sat, 11 Jul 2026 16:48:17 -0400 Subject: [PATCH 03/15] Open individual blocks and bypass cellimg caching layer --- .../spark/multisem/Wafer6061Inpainter.java | 108 +++++++++++------- 1 file changed, 66 insertions(+), 42 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 4452f30e5..99719bb6c 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -235,12 +235,36 @@ public void run() throws IOException { throw new IllegalArgumentException("expected a 3D tissue volume but " + params.fullDataset() + " has " + numDimensions + " dimensions"); } - if (n5.getDatasetAttributes(params.mask) == null) { + final DatasetAttributes maskAttributes = n5.getDatasetAttributes(params.mask); + if (maskAttributes == null) { throw new IllegalArgumentException("mask dataset " + params.mask + " does not exist"); } + // The inpainter reads and writes one chunk at a time (no whole-volume open, so no accumulating cell + // cache). That requires z to be a single chunk, so that every z-1 / z+1 section the z-average needs lives + // inside the block; and it reads the mask by the tissue block's grid position, so the two datasets must + // share the same block grid and origin. + if (s0Attributes.getBlockSize()[2] < s0Attributes.getDimensions()[2]) { + throw new IllegalArgumentException( + "block-local inpainting requires z to be a single chunk, but " + params.fullDataset() + + " has blockSize[2]=" + s0Attributes.getBlockSize()[2] + + " < dimensions[2]=" + s0Attributes.getDimensions()[2]); + } + if (! Arrays.equals(maskAttributes.getBlockSize(), s0Attributes.getBlockSize()) || + ! Arrays.equals(maskAttributes.getDimensions(), s0Attributes.getDimensions())) { + throw new IllegalArgumentException( + "mask grid " + Arrays.toString(maskAttributes.getDimensions()) + " @ " + + Arrays.toString(maskAttributes.getBlockSize()) + " must match tissue grid " + + Arrays.toString(s0Attributes.getDimensions()) + " @ " + + Arrays.toString(s0Attributes.getBlockSize())); + } + tissueTranslate = readTranslate(n5, params.fullDataset(), numDimensions); maskTranslate = readTranslate(n5, params.mask, numDimensions); + if (! Arrays.equals(tissueTranslate, maskTranslate)) { + throw new IllegalArgumentException("tissue translate " + Arrays.toString(tissueTranslate) + + " must match mask translate " + Arrays.toString(maskTranslate)); + } levels.add(params.fullDataset()); levelAttributes.put(params.fullDataset(), s0Attributes); @@ -272,15 +296,12 @@ public void run() throws IOException { final SparkConf conf = new SparkConf().setAppName(getClass().getSimpleName()); try (final JavaSparkContext sparkContext = new JavaSparkContext(conf)) { LOG.info("run: appId is {}", sparkContext.getConf().getAppId()); - runWithSparkContext(sparkContext, cloud, tissueTranslate, maskTranslate, - levels, levelAttributes, backupPath); + runWithSparkContext(sparkContext, cloud, levels, levelAttributes, backupPath); } } private void runWithSparkContext(final JavaSparkContext sparkContext, final PointCloud cloud, - final long[] tissueTranslate, - final long[] maskTranslate, final List levels, final Map levelAttributes, final String backupPath) { @@ -298,8 +319,6 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, blockIterator -> inpaintPartition(blockIterator, cloudBroadcast.value(), paramsBroadcast.value(), - tissueTranslate, - maskTranslate, backup)); final List modifiedS0 = modifiedRDD.collect(); @@ -552,8 +571,6 @@ private static RandomAccessibleInterval openDoubles(final N5Reader n private static Iterator inpaintPartition(final Iterator blocks, final PointCloud cloud, final Parameters params, - final long[] tissueTranslate, - final long[] maskTranslate, final String backupPath) { LogUtilities.setupExecutorLog4j("inpaint"); @@ -568,14 +585,7 @@ private static Iterator inpaintPartition(final Iterator bloc new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { final DatasetAttributes s0Attributes = reader.getDatasetAttributes(params.fullDataset()); - - final Img rawTissue = N5Utils.open(reader, params.fullDataset()); - final Img rawMask = N5Utils.open(reader, params.mask); - final RandomAccessibleInterval tissue = Views.translate(rawTissue, tissueTranslate); - final RandomAccessibleInterval mask = Views.translate(rawMask, maskTranslate); - - final long zMin = tissueTranslate[2]; - final long zMax = tissueTranslate[2] + s0Attributes.getDimensions()[2] - 1; + final DatasetAttributes maskAttributes = reader.getDatasetAttributes(params.mask); long considered = 0; long nearRoi = 0; @@ -595,22 +605,30 @@ private static Iterator inpaintPartition(final Iterator bloc } nearRoi++; - // (2) presence check: readBlock returns null for absent (empty) blocks. - final DataBlock originalBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); - if (originalBlock == null) { + // (2) presence check: readBlock returns null for absent (empty) blocks. The block it returns is also + // the raw tissue we inpaint from and back up, so no whole-volume open (and no accumulating cell cache) + // is needed: because z is a single chunk (guarded on the driver), every z-1 / z+1 section the + // z-average reads lives inside this same block. + final DataBlock tissueBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); + if (tissueBlock == null) { continue; } present++; - // (3) inpaint: fill mask==0 pixels with the z-average of the sections above/below. - final InpaintResult result = inpaintBlock(tissue, mask, block, tissueTranslate, zMin, zMax); + // (3) inpaint: fill mask==0 pixels with the z-average of the sections above/below. The mask is read by + // the same grid position (the datasets share a grid, guarded on the driver); an absent mask block + // counts as all-background (all holes), matching a zero-filled whole-volume read of a missing chunk. + final DataBlock maskBlock = reader.readBlock(params.mask, maskAttributes, block.gridPosition); + final InpaintResult result = inpaintBlock(asByteImg(tissueBlock, block.dimensions), + asByteImg(maskBlock, block.dimensions), + block.dimensions); if (! result.changed) { continue; } modified.add(block.gridPosition); if (! params.dryRun) { - backupWriter.writeBlock(params.fullDataset(), s0Attributes, originalBlock); + backupWriter.writeBlock(params.fullDataset(), s0Attributes, tissueBlock); N5Utils.saveBlock(result.inpainted, tissueWriter, params.fullDataset(), s0Attributes, block.gridPosition); } LOG.info("inpaintPartition: inpainted block {} (roiDistance={})", @@ -623,40 +641,46 @@ private static Iterator inpaintPartition(final Iterator bloc return modified.iterator(); } + /** + * Wraps a uint8 {@link DataBlock} as a block-local image. An absent (null) block becomes all-zeros, matching the + * zero fill an {@code N5Utils.open} whole-volume read would give for a missing chunk. + */ + private static Img asByteImg(final DataBlock dataBlock, final long[] blockDimensions) { + if (dataBlock == null) { + return ArrayImgs.unsignedBytes(blockDimensions); + } + return ArrayImgs.unsignedBytes((byte[]) dataBlock.getData(), blockDimensions); + } + /** * Produces the inpainted version of a single block: every pixel where the mask is 0 is replaced with the average - * of the tissue in the sections above and below (the single neighbor is copied at the volume's z-boundary). The - * {@code tissue} and {@code mask} are expected to be world-translated (i.e. indexed in the same coordinate frame, - * offset by their respective {@code translate}); {@code zMin}/{@code zMax} are the world z-bounds of the tissue. + * of the tissue in the sections above and below (the single neighbor is copied at the block's z-boundary). Both + * {@code tissueBlock} and {@code maskBlock} are block-local images with dimensions {@code blockDimensions}; because + * the volume is a single z-chunk, the block's z-boundary is the volume's z-boundary, so no neighbouring block is + * needed for the z-average. */ - static InpaintResult inpaintBlock(final RandomAccessibleInterval tissue, - final RandomAccessibleInterval mask, - final Grid.Block block, - final long[] tissueTranslate, - final long zMin, - final long zMax) { + static InpaintResult inpaintBlock(final RandomAccessibleInterval tissueBlock, + final RandomAccessibleInterval maskBlock, + final long[] blockDimensions) { - final RandomAccess tissueAccess = tissue.randomAccess(); - final RandomAccess maskAccess = Views.extendZero(mask).randomAccess(); + final RandomAccess tissueAccess = tissueBlock.randomAccess(); + final RandomAccess maskAccess = maskBlock.randomAccess(); + final long zMax = blockDimensions[2] - 1; - final Img inpainted = ArrayImgs.unsignedBytes(block.dimensions); + final Img inpainted = ArrayImgs.unsignedBytes(blockDimensions); final Cursor cursor = inpainted.localizingCursor(); final long[] local = new long[3]; - final long[] world = new long[3]; boolean changed = false; while (cursor.hasNext()) { final UnsignedByteType target = cursor.next(); cursor.localize(local); - world[0] = block.offset[0] + tissueTranslate[0] + local[0]; - world[1] = block.offset[1] + tissueTranslate[1] + local[1]; - world[2] = block.offset[2] + tissueTranslate[2] + local[2]; - final int original = tissueAccess.setPositionAndGet(world).get(); + final int original = tissueAccess.setPositionAndGet(local).get(); final int value; - if (maskAccess.setPositionAndGet(world).get() > 0) { + if (maskAccess.setPositionAndGet(local).get() > 0) { value = original; } else { - value = zAverage(tissueAccess, world, zMin, zMax); + value = zAverage(tissueAccess, local, 0, zMax); if (value != original) { changed = true; } From 05419a7f9cf2d588fb9762c042ef25fa642b3337 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sat, 11 Jul 2026 17:00:07 -0400 Subject: [PATCH 04/15] Add more diagnostics --- .../spark/multisem/Wafer6061Inpainter.java | 146 +++++++++++++++--- 1 file changed, 128 insertions(+), 18 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 99719bb6c..420b040e7 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -3,8 +3,11 @@ import com.beust.jcommander.Parameter; +import java.io.BufferedWriter; import java.io.IOException; import java.io.Serializable; +import java.nio.file.Files; +import java.nio.file.Paths; import java.util.ArrayList; import java.util.Arrays; import java.util.Iterator; @@ -147,6 +150,14 @@ public static class Parameters extends CommandLineParameters { description = "Only compute and log candidate / would-be-modified counts; do not write or back up anything.") public boolean dryRun = false; + @Parameter( + names = "--diagnosticsPath", + description = "Optional CSV file for per-block decisions (outside_roi / inside_roi / inpainted, plus the " + + "minimal inpainted z-layer). When set, decisions for ALL blocks are collected. Near-ROI blocks " + + "with no tissue are omitted (the distance filter runs before any presence check, so block " + + "emptiness is not observable outside the ROI and is not recorded as a category).") + public String diagnosticsPath; + public String fullDataset() { return dataset + "/s0"; } @@ -313,17 +324,37 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, final Broadcast cloudBroadcast = sparkContext.broadcast(cloud); final Broadcast paramsBroadcast = sparkContext.broadcast(params); - // 2. + 3. Filter (distance, then presence) and inpaint in one distributed pass. + // 2. + 3. Filter (distance, then presence) and inpaint in one distributed pass. When a diagnostics file is + // requested, every block's decision is collected (not just the inpainted ones) so the full picture can be drawn. final String backup = params.dryRun ? null : backupPath; - final JavaRDD modifiedRDD = sparkContext.parallelize(s0Blocks).mapPartitions( + final boolean emitAllDecisions = params.diagnosticsPath != null; + final JavaRDD decisionRDD = sparkContext.parallelize(s0Blocks).mapPartitions( blockIterator -> inpaintPartition(blockIterator, cloudBroadcast.value(), paramsBroadcast.value(), - backup)); + backup, + emitAllDecisions)); + + final List decisions = decisionRDD.collect(); - final List modifiedS0 = modifiedRDD.collect(); + final List modifiedS0 = new ArrayList<>(); + long outsideCount = 0; + long insideCount = 0; + for (final BlockDecision d : decisions) { + switch (d.decision) { + case DECISION_OUTSIDE_ROI: outsideCount++; break; + case DECISION_INSIDE_ROI: insideCount++; break; + case DECISION_INPAINTED: modifiedS0.add(new long[] {d.gridX, d.gridY, 0}); break; + default: break; + } + } LOG.info("runWithSparkContext: {} s0 block(s) were {}", modifiedS0.size(), params.dryRun ? "identified for inpainting (dry run)" : "inpainted"); + if (emitAllDecisions) { + LOG.info("runWithSparkContext: decision counts: outside_roi={}, inside_roi={}, inpainted={}", + outsideCount, insideCount, modifiedS0.size()); + writeDiagnostics(params.diagnosticsPath, decisions, s0Attributes); + } if (params.dryRun || modifiedS0.isEmpty()) { return; @@ -333,6 +364,32 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, updatePyramid(sparkContext, levels, levelAttributes, modifiedS0, backupPath); } + /** + * Writes one CSV row per recorded {@link BlockDecision} plus a small metadata header. Consumed as-is by the + * visualization; the client is the single source of truth for the decisions (no logic is re-derived downstream). + */ + private void writeDiagnostics(final String diagnosticsPath, + final List decisions, + final DatasetAttributes s0Attributes) { + final long[] dims = s0Attributes.getDimensions(); + final int[] blockSize = s0Attributes.getBlockSize(); + final long gridX = (dims[0] + blockSize[0] - 1) / blockSize[0]; + final long gridY = (dims[1] + blockSize[1] - 1) / blockSize[1]; + try (final BufferedWriter writer = Files.newBufferedWriter(Paths.get(diagnosticsPath))) { + writer.write("# Wafer6061Inpainter diagnostics\n"); + writer.write("# serial=" + params.serial + " maxRoiDistance=" + params.maxRoiDistance + + " gridX=" + gridX + " gridY=" + gridY + " zLayers=" + dims[2] + "\n"); + writer.write("gridX,gridY,roiDistance,decision,minLayer\n"); + for (final BlockDecision d : decisions) { + writer.write(d.gridX + "," + d.gridY + "," + Double.toString(d.roiDistance) + "," + + decisionLabel(d.decision) + "," + d.minLayer + "\n"); + } + } catch (final IOException e) { + throw new RuntimeException("failed to write diagnostics to " + diagnosticsPath, e); + } + LOG.info("writeDiagnostics: wrote {} block decisions to {}", decisions.size(), diagnosticsPath); + } + // ------------------------------------------------------------------------------------------------ // Step 1: load the per-slab 2-D point cloud from the xlog zarr. // ------------------------------------------------------------------------------------------------ @@ -568,14 +625,15 @@ private static RandomAccessibleInterval openDoubles(final N5Reader n // Steps 2 + 3: filter and inpaint one partition of s0 blocks. // ------------------------------------------------------------------------------------------------ - private static Iterator inpaintPartition(final Iterator blocks, - final PointCloud cloud, - final Parameters params, - final String backupPath) { + private static Iterator inpaintPartition(final Iterator blocks, + final PointCloud cloud, + final Parameters params, + final String backupPath, + final boolean emitAllDecisions) { LogUtilities.setupExecutorLog4j("inpaint"); - final List modified = new ArrayList<>(); + final List decisions = new ArrayList<>(); final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(params.k); try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); @@ -590,9 +648,12 @@ private static Iterator inpaintPartition(final Iterator bloc long considered = 0; long nearRoi = 0; long present = 0; + long inpainted = 0; while (blocks.hasNext()) { final Grid.Block block = blocks.next(); considered++; + final long gridX = block.gridPosition[0]; + final long gridY = block.gridPosition[1]; // (1) ROI-distance filter, computed with no I/O. The cloud is already in the s0 voxel frame, so the // block center (voxel index) queries it directly; distance_roi is interpolated by inverse-distance @@ -601,6 +662,9 @@ private static Iterator inpaintPartition(final Iterator bloc final double centerY = block.offset[1] + block.dimensions[1] / 2.0; final double roiDistance = cloud.interpolate(search, centerX, centerY, params.idwPower); if (! (roiDistance < params.maxRoiDistance)) { + if (emitAllDecisions) { + decisions.add(new BlockDecision(gridX, gridY, roiDistance, DECISION_OUTSIDE_ROI, -1)); + } continue; } nearRoi++; @@ -608,7 +672,7 @@ private static Iterator inpaintPartition(final Iterator bloc // (2) presence check: readBlock returns null for absent (empty) blocks. The block it returns is also // the raw tissue we inpaint from and back up, so no whole-volume open (and no accumulating cell cache) // is needed: because z is a single chunk (guarded on the driver), every z-1 / z+1 section the - // z-average reads lives inside this same block. + // z-average reads lives inside this same block. Near-ROI empty blocks are not recorded (see BlockDecision). final DataBlock tissueBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); if (tissueBlock == null) { continue; @@ -623,22 +687,26 @@ private static Iterator inpaintPartition(final Iterator bloc asByteImg(maskBlock, block.dimensions), block.dimensions); if (! result.changed) { + if (emitAllDecisions) { + decisions.add(new BlockDecision(gridX, gridY, roiDistance, DECISION_INSIDE_ROI, -1)); + } continue; } + inpainted++; - modified.add(block.gridPosition); + decisions.add(new BlockDecision(gridX, gridY, roiDistance, DECISION_INPAINTED, result.minChangedLayer)); if (! params.dryRun) { backupWriter.writeBlock(params.fullDataset(), s0Attributes, tissueBlock); N5Utils.saveBlock(result.inpainted, tissueWriter, params.fullDataset(), s0Attributes, block.gridPosition); } - LOG.info("inpaintPartition: inpainted block {} (roiDistance={})", - Arrays.toString(block.gridPosition), roiDistance); + LOG.info("inpaintPartition: inpainted block {} at min z-layer {} (roiDistance={})", + Arrays.toString(block.gridPosition), result.minChangedLayer, roiDistance); } - LOG.info("inpaintPartition: partition summary: considered={}, nearRoi(<{}um)={}, present={}, withHoles(modified)={}", - considered, params.maxRoiDistance, nearRoi, present, modified.size()); + LOG.info("inpaintPartition: partition summary: considered={}, nearRoi(<{}um)={}, present={}, inpainted={}", + considered, params.maxRoiDistance, nearRoi, present, inpainted); } - return modified.iterator(); + return decisions.iterator(); } /** @@ -671,6 +739,7 @@ static InpaintResult inpaintBlock(final RandomAccessibleInterval cursor = inpainted.localizingCursor(); final long[] local = new long[3]; boolean changed = false; + int minChangedLayer = Integer.MAX_VALUE; while (cursor.hasNext()) { final UnsignedByteType target = cursor.next(); cursor.localize(local); @@ -683,12 +752,15 @@ static InpaintResult inpaintBlock(final RandomAccessibleInterval inpainted; final boolean changed; + final int minChangedLayer; // minimal z-layer with an inpainted (changed) pixel, or -1 when nothing changed - InpaintResult(final Img inpainted, final boolean changed) { + InpaintResult(final Img inpainted, final boolean changed, final int minChangedLayer) { this.inpainted = inpainted; this.changed = changed; + this.minChangedLayer = minChangedLayer; + } + } + + // Per-block decision categories recorded for diagnostics / visualization. + static final int DECISION_OUTSIDE_ROI = 0; + static final int DECISION_INSIDE_ROI = 1; + static final int DECISION_INPAINTED = 2; + + static String decisionLabel(final int decision) { + switch (decision) { + case DECISION_OUTSIDE_ROI: return "outside_roi"; + case DECISION_INSIDE_ROI: return "inside_roi"; + case DECISION_INPAINTED: return "inpainted"; + default: return "unknown"; + } + } + + /** + * The recorded decision for a single s0 block (its z grid position is always 0 — z is a single chunk). Near-ROI + * blocks with no tissue are intentionally NOT represented: the distance filter runs before any presence check, so + * block emptiness is only observable inside the ROI and is not treated as a category. + */ + static class BlockDecision implements Serializable { + final long gridX; + final long gridY; + final double roiDistance; + final int decision; + final int minLayer; // minimal inpainted z-layer, or -1 when the block was not inpainted + + BlockDecision(final long gridX, final long gridY, final double roiDistance, + final int decision, final int minLayer) { + this.gridX = gridX; + this.gridY = gridY; + this.roiDistance = roiDistance; + this.decision = decision; + this.minLayer = minLayer; } } From 02fde8cdadc61a56a02934fa293d0969b1057c8e Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sat, 11 Jul 2026 17:08:01 -0400 Subject: [PATCH 05/15] Simplify logging behavior --- .../spark/multisem/Wafer6061Inpainter.java | 153 +++++------------- 1 file changed, 36 insertions(+), 117 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 420b040e7..d7516d108 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -3,11 +3,8 @@ import com.beust.jcommander.Parameter; -import java.io.BufferedWriter; import java.io.IOException; import java.io.Serializable; -import java.nio.file.Files; -import java.nio.file.Paths; import java.util.ArrayList; import java.util.Arrays; import java.util.Iterator; @@ -150,14 +147,6 @@ public static class Parameters extends CommandLineParameters { description = "Only compute and log candidate / would-be-modified counts; do not write or back up anything.") public boolean dryRun = false; - @Parameter( - names = "--diagnosticsPath", - description = "Optional CSV file for per-block decisions (outside_roi / inside_roi / inpainted, plus the " + - "minimal inpainted z-layer). When set, decisions for ALL blocks are collected. Near-ROI blocks " + - "with no tissue are omitted (the distance filter runs before any presence check, so block " + - "emptiness is not observable outside the ROI and is not recorded as a category).") - public String diagnosticsPath; - public String fullDataset() { return dataset + "/s0"; } @@ -324,37 +313,25 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, final Broadcast cloudBroadcast = sparkContext.broadcast(cloud); final Broadcast paramsBroadcast = sparkContext.broadcast(params); - // 2. + 3. Filter (distance, then presence) and inpaint in one distributed pass. When a diagnostics file is - // requested, every block's decision is collected (not just the inpainted ones) so the full picture can be drawn. + // A single driver line carries the grid extent and z-depth the visualization needs; the per-block + // 'blockDecision ...' lines are logged on the executors (see inpaintPartition / logDecision). + final long[] dims = s0Attributes.getDimensions(); + final int[] blockSize = s0Attributes.getBlockSize(); + LOG.info("runWithSparkContext: diagnostics metadata serial={} maxRoiDistance={} gridX={} gridY={} zLayers={}", + params.serial, params.maxRoiDistance, + (dims[0] + blockSize[0] - 1) / blockSize[0], (dims[1] + blockSize[1] - 1) / blockSize[1], dims[2]); + + // 2. + 3. Filter (distance, then presence) and inpaint in one distributed pass. final String backup = params.dryRun ? null : backupPath; - final boolean emitAllDecisions = params.diagnosticsPath != null; - final JavaRDD decisionRDD = sparkContext.parallelize(s0Blocks).mapPartitions( + final JavaRDD modifiedRDD = sparkContext.parallelize(s0Blocks).mapPartitions( blockIterator -> inpaintPartition(blockIterator, cloudBroadcast.value(), paramsBroadcast.value(), - backup, - emitAllDecisions)); - - final List decisions = decisionRDD.collect(); + backup)); - final List modifiedS0 = new ArrayList<>(); - long outsideCount = 0; - long insideCount = 0; - for (final BlockDecision d : decisions) { - switch (d.decision) { - case DECISION_OUTSIDE_ROI: outsideCount++; break; - case DECISION_INSIDE_ROI: insideCount++; break; - case DECISION_INPAINTED: modifiedS0.add(new long[] {d.gridX, d.gridY, 0}); break; - default: break; - } - } + final List modifiedS0 = modifiedRDD.collect(); LOG.info("runWithSparkContext: {} s0 block(s) were {}", modifiedS0.size(), params.dryRun ? "identified for inpainting (dry run)" : "inpainted"); - if (emitAllDecisions) { - LOG.info("runWithSparkContext: decision counts: outside_roi={}, inside_roi={}, inpainted={}", - outsideCount, insideCount, modifiedS0.size()); - writeDiagnostics(params.diagnosticsPath, decisions, s0Attributes); - } if (params.dryRun || modifiedS0.isEmpty()) { return; @@ -364,32 +341,6 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, updatePyramid(sparkContext, levels, levelAttributes, modifiedS0, backupPath); } - /** - * Writes one CSV row per recorded {@link BlockDecision} plus a small metadata header. Consumed as-is by the - * visualization; the client is the single source of truth for the decisions (no logic is re-derived downstream). - */ - private void writeDiagnostics(final String diagnosticsPath, - final List decisions, - final DatasetAttributes s0Attributes) { - final long[] dims = s0Attributes.getDimensions(); - final int[] blockSize = s0Attributes.getBlockSize(); - final long gridX = (dims[0] + blockSize[0] - 1) / blockSize[0]; - final long gridY = (dims[1] + blockSize[1] - 1) / blockSize[1]; - try (final BufferedWriter writer = Files.newBufferedWriter(Paths.get(diagnosticsPath))) { - writer.write("# Wafer6061Inpainter diagnostics\n"); - writer.write("# serial=" + params.serial + " maxRoiDistance=" + params.maxRoiDistance + - " gridX=" + gridX + " gridY=" + gridY + " zLayers=" + dims[2] + "\n"); - writer.write("gridX,gridY,roiDistance,decision,minLayer\n"); - for (final BlockDecision d : decisions) { - writer.write(d.gridX + "," + d.gridY + "," + Double.toString(d.roiDistance) + "," + - decisionLabel(d.decision) + "," + d.minLayer + "\n"); - } - } catch (final IOException e) { - throw new RuntimeException("failed to write diagnostics to " + diagnosticsPath, e); - } - LOG.info("writeDiagnostics: wrote {} block decisions to {}", decisions.size(), diagnosticsPath); - } - // ------------------------------------------------------------------------------------------------ // Step 1: load the per-slab 2-D point cloud from the xlog zarr. // ------------------------------------------------------------------------------------------------ @@ -625,15 +576,14 @@ private static RandomAccessibleInterval openDoubles(final N5Reader n // Steps 2 + 3: filter and inpaint one partition of s0 blocks. // ------------------------------------------------------------------------------------------------ - private static Iterator inpaintPartition(final Iterator blocks, - final PointCloud cloud, - final Parameters params, - final String backupPath, - final boolean emitAllDecisions) { + private static Iterator inpaintPartition(final Iterator blocks, + final PointCloud cloud, + final Parameters params, + final String backupPath) { LogUtilities.setupExecutorLog4j("inpaint"); - final List decisions = new ArrayList<>(); + final List modified = new ArrayList<>(); final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(params.k); try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); @@ -648,7 +598,6 @@ private static Iterator inpaintPartition(final Iterator inpaintPartition(final Iterator inpaintPartition(final Iterator tissueBlock = reader.readBlock(params.fullDataset(), s0Attributes, block.gridPosition); if (tissueBlock == null) { continue; @@ -687,26 +635,33 @@ private static Iterator inpaintPartition(final Iterator Date: Sun, 12 Jul 2026 15:35:33 -0400 Subject: [PATCH 06/15] Reduce the number of cli parameters --- .../spark/multisem/Wafer6061Inpainter.java | 172 ++++++++++-------- 1 file changed, 92 insertions(+), 80 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index d7516d108..fa3cc2fbb 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -12,6 +12,8 @@ import java.util.List; import java.util.Map; import java.util.function.Consumer; +import java.util.regex.Matcher; +import java.util.regex.Pattern; import bdv.export.Downsample; import net.imglib2.Cursor; @@ -68,6 +70,18 @@ */ public class Wafer6061Inpainter { + // SFOV image size fallback (pixels), used only when the xlog lacks x_sfov / y_sfov; the xlog values are + // preferred at runtime (see loadPointCloud), so these are wafer 60/61 defaults, not a tuning knob. + private static final int DEFAULT_SFOV_WIDTH = 2000; + private static final int DEFAULT_SFOV_HEIGHT = 1748; + + // Inverse-distance-weighting knobs for the ROI-distance interpolation; fixed for wafer 60/61. + private static final int IDW_K = 8; + private static final double IDW_POWER = 2.0; + + // Tissue containers are named w_s_r (e.g. w61_s109_r00); the serial is parsed from that. + private static final Pattern SERIAL_IN_NAME = Pattern.compile("_s(\\d+)"); + private static final Logger LOG = LoggerFactory.getLogger(Wafer6061Inpainter.class); public static class Parameters extends CommandLineParameters { @@ -97,46 +111,16 @@ public static class Parameters extends CommandLineParameters { @Parameter( names = "--serial", - description = "Serial label (id_serial) of the section stored in this N5; resolved to the reference arrays' slab position.", - required = true) - public long serial; - - @Parameter( - names = "--sfovWidth", - description = "SFOV image width in pixels (xlog X_SFOV size); used for the ingestion half-SFOV placement offset.") - public int sfovWidth = 2000; - - @Parameter( - names = "--sfovHeight", - description = "SFOV image height in pixels (xlog Y_SFOV size); used for the ingestion half-SFOV placement offset.") - public int sfovHeight = 1748; - - @Parameter( - names = "--scan", - description = "xlog scan index whose x/y positions define the cloud. Default (-1) auto-picks the first scan " + - "with finite x and rotation_slab for the slab (positions are nearly scan-independent).") - public int scan = -1; + description = "Serial label (id_serial) of the section stored in this N5, resolved to the reference " + + "arrays' slab position. Defaults to the serial parsed from the container name " + + "(w_s_r, e.g. w61_s109_r00 -> 109).") + public long serial = -1; @Parameter( names = "--maxRoiDistance", description = "Keep (inpaint) a block only if its interpolated distance_roi (microns) is less than this.") public double maxRoiDistance = 10.0; - @Parameter( - names = "--k", - description = "Number of nearest neighbors used for inverse-distance weighting.") - public int k = 8; - - @Parameter( - names = "--idwPower", - description = "Power p in the inverse-distance weights 1 / r^p.") - public double idwPower = 2.0; - - @Parameter( - names = "--downsampleFactors", - description = "Relative per-level downsampling factors of the existing pyramid, e.g. 2,2,1.") - public String downsampleFactors = "2,2,1"; - @Parameter( names = "--backupPath", description = "N5 container where original (overwritten) blocks are backed up. Defaults to a sibling _backup.n5.") @@ -151,15 +135,6 @@ public String fullDataset() { return dataset + "/s0"; } - public int[] getDownsampleFactors() { - final String[] parts = downsampleFactors.split(","); - final int[] factors = new int[parts.length]; - for (int i = 0; i < parts.length; i++) { - factors[i] = Integer.parseInt(parts[i].trim()); - } - return factors; - } - public String getBackupPath() { if (backupPath != null) { return backupPath; @@ -175,16 +150,31 @@ public String getBackupPath() { } public void validate() { - if (sfovWidth <= 0 || sfovHeight <= 0) { - throw new IllegalArgumentException("--sfovWidth and --sfovHeight must be positive"); - } if (maxRoiDistance <= 0) { throw new IllegalArgumentException("--maxRoiDistance must be positive"); } - if (k < 1) { - throw new IllegalArgumentException("--k must be at least 1"); + if (serial < 0) { + serial = inferSerial(n5Path); } } + + /** Parses the serial label from a container name like {@code .../w61_s109_r00} (basename {@code _s}). */ + static long inferSerial(final String n5Path) { + String basename = n5Path; + while (basename.endsWith("/")) { + basename = basename.substring(0, basename.length() - 1); + } + final int slash = basename.lastIndexOf('/'); + if (slash >= 0) { + basename = basename.substring(slash + 1); + } + final Matcher matcher = SERIAL_IN_NAME.matcher(basename); + if (matcher.find()) { + return Long.parseLong(matcher.group(1)); + } + throw new IllegalArgumentException("could not infer the serial from n5Path '" + n5Path + + "'; pass --serial explicitly"); + } } private final Parameters params; @@ -214,8 +204,7 @@ public void runClient(final String[] args) throws Exception { public void run() throws IOException { // 1. Load the (small) per-slab 2-D point cloud from the xlog on the driver, placed in the tissue voxel frame. - final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial, - params.scan, params.sfovWidth, params.sfovHeight); + final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial); LOG.info("run: loaded {} ROI reference points for serial {}", cloud.size(), params.serial); // Read the tissue s0 metadata and discover the existing pyramid levels. @@ -224,6 +213,7 @@ public void run() throws IOException { final int numDimensions; final List levels = new ArrayList<>(); final Map levelAttributes = new LinkedHashMap<>(); + int[] downsampleFactors = null; try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path)) { final DatasetAttributes s0Attributes = n5.getDatasetAttributes(params.fullDataset()); @@ -276,9 +266,17 @@ public void run() throws IOException { levels.add(levelDataset); levelAttributes.put(levelDataset, n5.getDatasetAttributes(levelDataset)); } + + // The relative per-level downsampling factor is read from the pyramid itself rather than passed in: s1's + // "downsamplingFactors" attribute is the factor relative to s0, and these pyramids are built with a + // constant factor at every step (see DownsampleHelper / N5DownsamplerSpark), so it applies to all levels. + if (levels.size() > 1) { + downsampleFactors = readDownsamplingFactors(n5, levels.get(1), numDimensions); + } } - LOG.info("run: tissue translate={}, mask translate={}, pyramid levels={}", - Arrays.toString(tissueTranslate), Arrays.toString(maskTranslate), levels); + LOG.info("run: tissue translate={}, mask translate={}, pyramid levels={}, downsampleFactors={}", + Arrays.toString(tissueTranslate), Arrays.toString(maskTranslate), levels, + Arrays.toString(downsampleFactors)); // Create the backup container and mirror all datasets that might receive backups. final String backupPath = params.getBackupPath(); @@ -296,7 +294,7 @@ public void run() throws IOException { final SparkConf conf = new SparkConf().setAppName(getClass().getSimpleName()); try (final JavaSparkContext sparkContext = new JavaSparkContext(conf)) { LOG.info("run: appId is {}", sparkContext.getConf().getAppId()); - runWithSparkContext(sparkContext, cloud, levels, levelAttributes, backupPath); + runWithSparkContext(sparkContext, cloud, levels, levelAttributes, backupPath, downsampleFactors); } } @@ -304,7 +302,8 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, final PointCloud cloud, final List levels, final Map levelAttributes, - final String backupPath) { + final String backupPath, + final int[] downsampleFactors) { final DatasetAttributes s0Attributes = levelAttributes.get(params.fullDataset()); final List s0Blocks = Grid.create(s0Attributes.getDimensions(), s0Attributes.getBlockSize()); @@ -338,7 +337,7 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, } // 6. Selectively update the downsample pyramid, one level at a time. - updatePyramid(sparkContext, levels, levelAttributes, modifiedS0, backupPath); + updatePyramid(sparkContext, levels, levelAttributes, modifiedS0, backupPath, downsampleFactors); } // ------------------------------------------------------------------------------------------------ @@ -360,10 +359,7 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, * positions slightly, well within {@code maxRoiDistance}, so the unaligned ingestion placement is used directly. */ static PointCloud loadPointCloud(final String xlogPath, - final long serial, - final int scanOverride, - final int sfovWidth, - final int sfovHeight) { + final long serial) { try (final N5Reader xlog = new N5Factory().openReader(xlogPath)) { final double[] idSerial = read1d(xlog, "id_serial"); @@ -378,9 +374,9 @@ static PointCloud loadPointCloud(final String xlogPath, } final long slabCount = idSerial.length; // 413 for wafer 61 - // SFOV image size defines the half-SFOV placement offset; prefer the xlog (x_sfov / y_sfov), fall back to args. - final int sfW = sfovSize(xlog, "x_sfov", sfovWidth); - final int sfH = sfovSize(xlog, "y_sfov", sfovHeight); + // SFOV image size defines the half-SFOV placement offset; prefer the xlog (x_sfov / y_sfov), fall back to constants. + final int sfW = sfovSize(xlog, "x_sfov", DEFAULT_SFOV_WIDTH); + final int sfH = sfovSize(xlog, "y_sfov", DEFAULT_SFOV_HEIGHT); // distance_roi is scan-independent: [slab, mfov, sfov] -> 2-D (sfov, mfov) after slicing the slab axis. final RandomAccessibleInterval distSlab = readSlab(xlog, "distance_roi", slabCount, slabPosition); @@ -395,23 +391,22 @@ static PointCloud loadPointCloud(final String xlogPath, final int slabAxisX = axisOfSize(xAll, slabCount, "x"); final int scanAxisX = axisOfSize(xAll, nScans, "x"); - // Choose the scan whose x/y define the cloud (positions are nearly scan-independent). - int scan = scanOverride; - if (scan < 0) { - for (int s = 0; s < nScans; s++) { - if (! Double.isFinite(rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, s))) { - continue; - } - if (hasFinite(sliceScanAndSlab(xAll, scanAxisX, s, slabAxisX, slabPosition))) { - scan = s; - break; - } + // Choose the scan whose x/y define the cloud: the first with finite x and rotation_slab for this slab. + // The positions are nearly scan-independent, so no scan override is needed. + int scan = -1; + for (int s = 0; s < nScans; s++) { + if (! Double.isFinite(rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, s))) { + continue; } - if (scan < 0) { - throw new IllegalArgumentException("no scan with finite x and rotation_slab found for serial " + - serial + " (slab position " + slabPosition + ")"); + if (hasFinite(sliceScanAndSlab(xAll, scanAxisX, s, slabAxisX, slabPosition))) { + scan = s; + break; } } + if (scan < 0) { + throw new IllegalArgumentException("no scan with finite x and rotation_slab found for serial " + + serial + " (slab position " + slabPosition + ")"); + } final double rotationSlab = rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, scan); final double theta = Math.toRadians(180.0 + rotationSlab); @@ -584,7 +579,7 @@ private static Iterator inpaintPartition(final Iterator bloc LogUtilities.setupExecutorLog4j("inpaint"); final List modified = new ArrayList<>(); - final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(params.k); + final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(IDW_K); try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); final N5Writer tissueWriter = params.dryRun ? null : @@ -609,7 +604,7 @@ private static Iterator inpaintPartition(final Iterator bloc // weighting (the interpolated value is in microns regardless of the voxel-space query units). final double centerX = block.offset[0] + block.dimensions[0] / 2.0; final double centerY = block.offset[1] + block.dimensions[1] / 2.0; - final double roiDistance = cloud.interpolate(search, centerX, centerY, params.idwPower); + final double roiDistance = cloud.interpolate(search, centerX, centerY, IDW_POWER); if (! (roiDistance < params.maxRoiDistance)) { logDecision(gridX, gridY, "outside_roi", -1, roiDistance); continue; @@ -756,9 +751,9 @@ private void updatePyramid(final JavaSparkContext sparkContext, final List levels, final Map levelAttributes, final List modifiedS0, - final String backupPath) { + final String backupPath, + final int[] factors) { - final int[] factors = params.getDownsampleFactors(); List modifiedPrevious = modifiedS0; for (int scale = 1; scale < levels.size(); scale++) { @@ -866,6 +861,23 @@ private static long[] readTranslate(final N5Reader n5, final String dataset, fin return translate != null ? translate : new long[numDimensions]; } + /** + * Reads the {@code downsamplingFactors} attribute of a pyramid level (the factor relative to s0, as written by + * {@code N5DownsamplerSpark} / render's export). Fails fast when it is absent so the pyramid factor is never guessed. + */ + private static int[] readDownsamplingFactors(final N5Reader n5, final String dataset, final int numDimensions) { + final int[] factors = n5.getAttribute(dataset, "downsamplingFactors", int[].class); + if (factors == null) { + throw new IllegalArgumentException("dataset " + dataset + " has no 'downsamplingFactors' attribute to derive " + + "the pyramid downsampling factor from"); + } + if (factors.length != numDimensions) { + throw new IllegalArgumentException("dataset " + dataset + " has downsamplingFactors " + + Arrays.toString(factors) + " but the volume is " + numDimensions + "-dimensional"); + } + return factors; + } + /** Result of inpainting one block: the (block-sized) inpainted image and whether any pixel changed. */ static class InpaintResult { final Img inpainted; From 14f228a560dcd12c340b5ef3de5c206e57b048e0 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sun, 12 Jul 2026 15:39:57 -0400 Subject: [PATCH 07/15] Simplify some other logic --- .../spark/multisem/Wafer6061Inpainter.java | 72 ++++++++----------- 1 file changed, 31 insertions(+), 41 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index fa3cc2fbb..9e07711e4 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -214,7 +214,7 @@ public void run() throws IOException { final List levels = new ArrayList<>(); final Map levelAttributes = new LinkedHashMap<>(); int[] downsampleFactors = null; - try (final N5Reader n5 = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path)) { + try (final N5Reader n5 = openN5Reader(params.n5Path)) { final DatasetAttributes s0Attributes = n5.getDatasetAttributes(params.fullDataset()); if (s0Attributes == null) { @@ -281,7 +281,7 @@ public void run() throws IOException { // Create the backup container and mirror all datasets that might receive backups. final String backupPath = params.getBackupPath(); if (! params.dryRun) { - try (final N5Writer backup = new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + try (final N5Writer backup = openN5Writer(backupPath)) { for (final String levelDataset : levels) { if (! backup.datasetExists(levelDataset)) { backup.createDataset(levelDataset, levelAttributes.get(levelDataset)); @@ -542,18 +542,7 @@ private static RandomAccessibleInterval readSlab(final N5Reader n5, final long slabCount, final int slabPosition) { final RandomAccessibleInterval img = openDoubles(n5, dataset); - int slabAxis = -1; - for (int d = 0; d < img.numDimensions(); d++) { - if (img.dimension(d) == slabCount) { - slabAxis = d; - break; - } - } - if (slabAxis < 0) { - throw new IllegalArgumentException("could not find slab axis (size " + slabCount + ") in " + dataset + - " with dimensions " + Arrays.toString(img.dimensionsAsLongArray())); - } - return Views.hyperSlice(img, slabAxis, slabPosition); + return Views.hyperSlice(img, axisOfSize(img, slabCount, dataset), slabPosition); } /** @@ -581,11 +570,9 @@ private static Iterator inpaintPartition(final Iterator bloc final List modified = new ArrayList<>(); final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(IDW_K); - try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, params.n5Path); - final N5Writer tissueWriter = params.dryRun ? null : - new N5Factory().openWriter(N5Factory.StorageFormat.N5, params.n5Path); - final N5Writer backupWriter = params.dryRun ? null : - new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + try (final N5Reader reader = openN5Reader(params.n5Path); + final N5Writer tissueWriter = params.dryRun ? null : openN5Writer(params.n5Path); + final N5Writer backupWriter = params.dryRun ? null : openN5Writer(backupPath)) { final DatasetAttributes s0Attributes = reader.getDatasetAttributes(params.fullDataset()); final DatasetAttributes maskAttributes = reader.getDatasetAttributes(params.mask); @@ -719,28 +706,21 @@ private static int zAverage(final RandomAccess tissueAccess, final long zMin, final long zMax) { final long z = world[2]; - final boolean hasAbove = (z - 1) >= zMin; - final boolean hasBelow = (z + 1) <= zMax; - if (hasAbove && hasBelow) { - world[2] = z - 1; - final int above = tissueAccess.setPositionAndGet(world).get(); - world[2] = z + 1; - final int below = tissueAccess.setPositionAndGet(world).get(); - world[2] = z; - return (above + below) >>> 1; - } else if (hasAbove) { + int sum = 0; + int count = 0; + if (z - 1 >= zMin) { world[2] = z - 1; - final int above = tissueAccess.setPositionAndGet(world).get(); - world[2] = z; - return above; - } else if (hasBelow) { + sum += tissueAccess.setPositionAndGet(world).get(); + count++; + } + if (z + 1 <= zMax) { world[2] = z + 1; - final int below = tissueAccess.setPositionAndGet(world).get(); - world[2] = z; - return below; - } else { - return tissueAccess.setPositionAndGet(world).get(); + sum += tissueAccess.setPositionAndGet(world).get(); + count++; } + world[2] = z; + // both neighbors -> mean (== the old (above+below)>>>1 for byte values); one -> that neighbor; none -> unchanged. + return count > 0 ? sum / count : tissueAccess.setPositionAndGet(world).get(); } // ------------------------------------------------------------------------------------------------ @@ -810,9 +790,9 @@ static void downsampleBlock(final long[] gridPosition, targetSize[d] = cellDimensions[d]; } - try (final N5Reader reader = new N5Factory().openReader(N5Factory.StorageFormat.N5, n5Path); - final N5Writer writer = new N5Factory().openWriter(N5Factory.StorageFormat.N5, n5Path); - final N5Writer backupWriter = new N5Factory().openWriter(N5Factory.StorageFormat.N5, backupPath)) { + try (final N5Reader reader = openN5Reader(n5Path); + final N5Writer writer = openN5Writer(n5Path); + final N5Writer backupWriter = openN5Writer(backupPath)) { final RandomAccessibleInterval source = N5Utils.open(reader, fromDataset); final RandomAccessibleInterval sourceBlock = Views.offsetInterval(source, sourceMin, sourceSize); @@ -861,6 +841,16 @@ private static long[] readTranslate(final N5Reader n5, final String dataset, fin return translate != null ? translate : new long[numDimensions]; } + /** Opens an N5 reader for a tissue/backup container (explicit N5 format; local path or gs://). */ + private static N5Reader openN5Reader(final String path) { + return new N5Factory().openReader(N5Factory.StorageFormat.N5, path); + } + + /** Opens an N5 writer for a tissue/backup container (explicit N5 format; local path or gs://). */ + private static N5Writer openN5Writer(final String path) { + return new N5Factory().openWriter(N5Factory.StorageFormat.N5, path); + } + /** * Reads the {@code downsamplingFactors} attribute of a pyramid level (the factor relative to s0, as written by * {@code N5DownsamplerSpark} / render's export). Fails fast when it is absent so the pyramid factor is never guessed. From 4a1d6d18d29fbd242f546375aaf516f67023f9bc Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sun, 12 Jul 2026 15:55:03 -0400 Subject: [PATCH 08/15] Hardcode xlog axes to make selection deterministic --- .../spark/multisem/Wafer6061Inpainter.java | 98 ++++++++++--------- 1 file changed, 50 insertions(+), 48 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 9e07711e4..52044c531 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -82,6 +82,23 @@ public class Wafer6061Inpainter { // Tissue containers are named w_s_r (e.g. w61_s109_r00); the serial is parsed from that. private static final Pattern SERIAL_IN_NAME = Pattern.compile("_s(\\d+)"); + // xlog field layout (consumed by loadPointCloud). The xlog is a zarr: its arrays are stored C-order, but the + // n5-zarr reader REVERSES the axis order, so the imglib2 axis indices below are the reverse of the C-order shape. + // field C-order shape imglib2 axes (what this code sees) meaning + // id_serial (slab,) [slab] serial label per slab position + // x, y (scan, slab, mfov, sfov) [sfov, mfov, slab, scan] SFOV centre, full-res px + // rotation_slab (scan, slab) [slab, scan] per-slab rotation, degrees + // distance_roi (slab, mfov, sfov) [sfov, mfov, slab] distance to ROI, um (scan-independent) + // x_sfov, y_sfov (width,), (height,) [size] SFOV image size, px + // The slab/scan axis indices are HARDCODED (below) rather than discovered by matching sizes: the sizes are not + // guaranteed to be unique. x / y / distance_roi share (mfov, sfov) indexing, so after slicing the slab (and scan) + // axis they all reduce to the same 2-D (sfov, mfov) plane with matching indices. + private static final int XY_AXIS_SLAB = 2; + private static final int XY_AXIS_SCAN = 3; + private static final int ROT_AXIS_SLAB = 0; + private static final int ROT_AXIS_SCAN = 1; + private static final int DIST_AXIS_SLAB = 2; + private static final Logger LOG = LoggerFactory.getLogger(Wafer6061Inpainter.class); public static class Parameters extends CommandLineParameters { @@ -378,27 +395,29 @@ static PointCloud loadPointCloud(final String xlogPath, final int sfW = sfovSize(xlog, "x_sfov", DEFAULT_SFOV_WIDTH); final int sfH = sfovSize(xlog, "y_sfov", DEFAULT_SFOV_HEIGHT); - // distance_roi is scan-independent: [slab, mfov, sfov] -> 2-D (sfov, mfov) after slicing the slab axis. - final RandomAccessibleInterval distSlab = readSlab(xlog, "distance_roi", slabCount, slabPosition); - - // x / y are [scan, slab, mfov, sfov]; rotation_slab is [scan, slab]. Identify axes by their (unique) sizes. + // Open the reference arrays (axis layout documented and hardcoded at the top of the class) and verify each + // one carries the slab count on its expected axis before we slice by it. final RandomAccessibleInterval xAll = openDoubles(xlog, "x"); final RandomAccessibleInterval yAll = openDoubles(xlog, "y"); final RandomAccessibleInterval rotAll = openDoubles(xlog, "rotation_slab"); - final int rotSlabAxis = axisOfSize(rotAll, slabCount, "rotation_slab"); - final int rotScanAxis = 1 - rotSlabAxis; - final long nScans = rotAll.dimension(rotScanAxis); - final int slabAxisX = axisOfSize(xAll, slabCount, "x"); - final int scanAxisX = axisOfSize(xAll, nScans, "x"); + final RandomAccessibleInterval distAll = openDoubles(xlog, "distance_roi"); + requireSlabAxis(xAll, XY_AXIS_SLAB, slabCount, "x"); + requireSlabAxis(yAll, XY_AXIS_SLAB, slabCount, "y"); + requireSlabAxis(rotAll, ROT_AXIS_SLAB, slabCount, "rotation_slab"); + requireSlabAxis(distAll, DIST_AXIS_SLAB, slabCount, "distance_roi"); + final long nScans = rotAll.dimension(ROT_AXIS_SCAN); + + // distance_roi is scan-independent: slice the slab axis -> 2-D (sfov, mfov). + final RandomAccessibleInterval distSlab = Views.hyperSlice(distAll, DIST_AXIS_SLAB, slabPosition); // Choose the scan whose x/y define the cloud: the first with finite x and rotation_slab for this slab. // The positions are nearly scan-independent, so no scan override is needed. int scan = -1; for (int s = 0; s < nScans; s++) { - if (! Double.isFinite(rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, s))) { + if (! Double.isFinite(rotationAt(rotAll, slabPosition, s))) { continue; } - if (hasFinite(sliceScanAndSlab(xAll, scanAxisX, s, slabAxisX, slabPosition))) { + if (hasFinite(sliceScanAndSlab(xAll, s, slabPosition))) { scan = s; break; } @@ -408,13 +427,13 @@ static PointCloud loadPointCloud(final String xlogPath, serial + " (slab position " + slabPosition + ")"); } - final double rotationSlab = rotationAt(rotAll, rotSlabAxis, rotScanAxis, slabPosition, scan); + final double rotationSlab = rotationAt(rotAll, slabPosition, scan); final double theta = Math.toRadians(180.0 + rotationSlab); final double cos = Math.cos(theta); final double sin = Math.sin(theta); - final RandomAccessibleInterval xSlab = sliceScanAndSlab(xAll, scanAxisX, scan, slabAxisX, slabPosition); - final RandomAccessibleInterval ySlab = sliceScanAndSlab(yAll, scanAxisX, scan, slabAxisX, slabPosition); + final RandomAccessibleInterval xSlab = sliceScanAndSlab(xAll, scan, slabPosition); + final RandomAccessibleInterval ySlab = sliceScanAndSlab(yAll, scan, slabPosition); // Place each SFOV center (rotation only) and keep the finite ones with their distance_roi. final List cxs = new ArrayList<>(); @@ -477,41 +496,33 @@ private static int sfovSize(final N5Reader xlog, final String dataset, final int return fallback; } - /** Index of the (unique-sized) axis matching {@code size}. */ - private static int axisOfSize(final RandomAccessibleInterval img, final long size, final String name) { - for (int d = 0; d < img.numDimensions(); d++) { - if (img.dimension(d) == size) { - return d; - } + /** Fails fast if a hardcoded slab axis does not carry the expected slab count (guards the layout at the top). */ + private static void requireSlabAxis(final RandomAccessibleInterval img, final int axis, final long slabCount, + final String field) { + if (img.numDimensions() <= axis || img.dimension(axis) != slabCount) { + throw new IllegalArgumentException("xlog field '" + field + "' has dimensions " + + Arrays.toString(img.dimensionsAsLongArray()) + " but the hardcoded layout " + + "expects " + slabCount + " slabs on axis " + axis + + " (see the xlog axis layout at the top of the class)"); } - throw new IllegalArgumentException("no axis of size " + size + " in " + name + " with dimensions " + - Arrays.toString(img.dimensionsAsLongArray())); } - /** Value of the 2-D {@code rotation_slab} at the given slab position and scan. */ + /** Value of the 2-D {@code rotation_slab} (axes {@code [slab, scan]}) at the given slab position and scan. */ private static double rotationAt(final RandomAccessibleInterval rotAll, - final int slabAxis, - final int scanAxis, - final int slabPosition, - final int scan) { + final long slabPosition, + final long scan) { final long[] p = new long[2]; - p[slabAxis] = slabPosition; - p[scanAxis] = scan; + p[ROT_AXIS_SLAB] = slabPosition; + p[ROT_AXIS_SCAN] = scan; return rotAll.randomAccess().setPositionAndGet(p).get(); } - /** Slices the 4-D {@code x}/{@code y} array to the 2-D (sfov, mfov) plane for one scan and slab. */ - private static RandomAccessibleInterval sliceScanAndSlab(final RandomAccessibleInterval img, - final int scanAxis, + /** Slices the 4-D {@code x}/{@code y} array (axes {@code [sfov, mfov, slab, scan]}) to the 2-D (sfov, mfov) plane. */ + private static RandomAccessibleInterval sliceScanAndSlab(final RandomAccessibleInterval xy, final long scan, - final int slabAxis, final long slabPosition) { - // hyper-slice the higher-indexed axis first so the lower index stays valid afterwards - final int hi = Math.max(scanAxis, slabAxis); - final int lo = Math.min(scanAxis, slabAxis); - final long hiPos = (hi == scanAxis) ? scan : slabPosition; - final long loPos = (lo == scanAxis) ? scan : slabPosition; - return Views.hyperSlice(Views.hyperSlice(img, hi, hiPos), lo, loPos); + // slice the higher axis (scan) first so the lower slab axis index stays valid. + return Views.hyperSlice(Views.hyperSlice(xy, XY_AXIS_SCAN, scan), XY_AXIS_SLAB, slabPosition); } /** True if the interval has at least one finite value. */ @@ -536,15 +547,6 @@ private static double[] read1d(final N5Reader n5, final String dataset) { return values; } - /** Opens a 3-D reference array and returns the 2-D slice for the given slab position. */ - private static RandomAccessibleInterval readSlab(final N5Reader n5, - final String dataset, - final long slabCount, - final int slabPosition) { - final RandomAccessibleInterval img = openDoubles(n5, dataset); - return Views.hyperSlice(img, axisOfSize(img, slabCount, dataset), slabPosition); - } - /** * Opens a double-typed dataset using the {@code (blockNotFoundHandler, accessFlags)} overload. This avoids * {@code N5Utils.open(reader, dataset)}'s {@code isLabelMultisetType -> getAttribute} path, which NPEs on blosc From a2238b79445bb3477f89b08e53ef13b3db31d51e Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sun, 12 Jul 2026 16:28:00 -0400 Subject: [PATCH 09/15] Load balance the inpainting pass better by shuffling --- .../spark/multisem/Wafer6061Inpainter.java | 110 ++++++++++++------ 1 file changed, 75 insertions(+), 35 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 52044c531..c2461a69f 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -7,10 +7,12 @@ import java.io.Serializable; import java.util.ArrayList; import java.util.Arrays; +import java.util.Collections; import java.util.Iterator; import java.util.LinkedHashMap; import java.util.List; import java.util.Map; +import java.util.Random; import java.util.function.Consumer; import java.util.regex.Matcher; import java.util.regex.Pattern; @@ -323,8 +325,18 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, final int[] downsampleFactors) { final DatasetAttributes s0Attributes = levelAttributes.get(params.fullDataset()); - final List s0Blocks = Grid.create(s0Attributes.getDimensions(), s0Attributes.getBlockSize()); - LOG.info("runWithSparkContext: {} s0 grid blocks to consider", s0Blocks.size()); + final List s0Blocks = new ArrayList<>(Grid.create(s0Attributes.getDimensions(), + s0Attributes.getBlockSize())); + + // Grid.create returns blocks in raster order, and Spark's parallelize slices the list into contiguous + // partitions. The blocks that actually do work (present-check + inpaint) are the near-ROI ones, and the ROI is + // a small, spatially clustered region, so contiguous slicing piles all the expensive blocks into a few + // partitions while the rest only run the cheap no-I/O distance filter -> severe load skew. Shuffling first gives + // every partition a representative mix of near- and far-ROI blocks, so the pass is balanced. The seed (serial) + // keeps the partitioning reproducible, and the outputs (emitted grid positions, per-block decision logs) are + // order-independent, so this changes only the distribution of work, not the result. + Collections.shuffle(s0Blocks, new Random(params.serial)); + LOG.info("runWithSparkContext: {} s0 grid blocks to consider (shuffled for load balance)", s0Blocks.size()); final Broadcast cloudBroadcast = sparkContext.broadcast(cloud); final Broadcast paramsBroadcast = sparkContext.broadcast(params); @@ -750,35 +762,70 @@ private void updatePyramid(final JavaSparkContext sparkContext, affected.putIfAbsent(Arrays.toString(g), g); } final List affectedBlocks = new ArrayList<>(affected.values()); + // Every affected block does real work, so there is no near/far skew here, but the per-block cost still varies + // with how much of its source region is actually present (dense in the ROI interior, sparse at its edge). + // The affected blocks inherit a spatial order, so contiguous partitioning would group same-density + // neighbours together and leave some partitions all-dense and others all-sparse. Shuffle (as for s0) mixes + // densities across partitions; the seed varies per level but stays reproducible. + Collections.shuffle(affectedBlocks, new Random(params.serial * 31L + scale)); LOG.info("updatePyramid: re-downsampling {} block(s) for {}", affectedBlocks.size(), toDataset); final String n5Path = params.n5Path; - sparkContext.parallelize(affectedBlocks).foreach( - gridPosition -> downsampleBlock(gridPosition, n5Path, backupPath, - fromDataset, toDataset, toAttributes, factors)); + sparkContext.parallelize(affectedBlocks).foreachPartition( + gridPositions -> downsamplePartition(gridPositions, n5Path, backupPath, + fromDataset, toDataset, toAttributes, factors)); modifiedPrevious = affectedBlocks; } } /** - * Re-downsamples a single pyramid block from the (already updated) previous level, replicating the per-block - * math of {@code N5DownsamplerSpark} so the result matches the rest of the pyramid. The original block is backed - * up before being overwritten. + * Re-downsamples one partition's worth of pyramid blocks, opening the N5 handles and the source level once + * for the whole partition rather than per block (the source open is a lazy {@code N5Utils.open}, so each block + * still reads only the source chunks it needs). Each block is rebuilt from the (already updated) previous level + * with the same per-block math as {@code N5DownsamplerSpark} so the result matches the rest of the pyramid, backing + * up the original block before overwriting it. */ - static void downsampleBlock(final long[] gridPosition, - final String n5Path, - final String backupPath, - final String fromDataset, - final String toDataset, - final DatasetAttributes toAttributes, - final int[] factors) { + static void downsamplePartition(final Iterator gridPositions, + final String n5Path, + final String backupPath, + final String fromDataset, + final String toDataset, + final DatasetAttributes toAttributes, + final int[] factors) { LogUtilities.setupExecutorLog4j("downsample"); - final int n = toAttributes.getNumDimensions(); final CellGrid cellGrid = new CellGrid(toAttributes.getDimensions(), toAttributes.getBlockSize()); + try (final N5Reader reader = openN5Reader(n5Path); + final N5Writer writer = openN5Writer(n5Path); + final N5Writer backupWriter = openN5Writer(backupPath)) { + + final RandomAccessibleInterval source = N5Utils.open(reader, fromDataset); + + int count = 0; + while (gridPositions.hasNext()) { + downsampleBlock(source, reader, writer, backupWriter, cellGrid, + gridPositions.next(), toDataset, toAttributes, factors); + count++; + } + LOG.info("downsamplePartition: re-downsampled {} {} block(s)", count, toDataset); + } + } + + /** Re-downsamples a single block using the handles and source level already opened for the partition. */ + private static void downsampleBlock(final RandomAccessibleInterval source, + final N5Reader reader, + final N5Writer writer, + final N5Writer backupWriter, + final CellGrid cellGrid, + final long[] gridPosition, + final String toDataset, + final DatasetAttributes toAttributes, + final int[] factors) { + + final int n = toAttributes.getNumDimensions(); final long[] targetMin = new long[n]; final int[] cellDimensions = new int[n]; cellGrid.getCellDimensions(gridPosition, targetMin, cellDimensions); @@ -792,27 +839,20 @@ static void downsampleBlock(final long[] gridPosition, targetSize[d] = cellDimensions[d]; } - try (final N5Reader reader = openN5Reader(n5Path); - final N5Writer writer = openN5Writer(n5Path); - final N5Writer backupWriter = openN5Writer(backupPath)) { - - final RandomAccessibleInterval source = N5Utils.open(reader, fromDataset); - final RandomAccessibleInterval sourceBlock = Views.offsetInterval(source, sourceMin, sourceSize); - - final Img targetBlock = ArrayImgs.unsignedBytes(targetSize); - Downsample.downsample(sourceBlock, targetBlock, factors); + final RandomAccessibleInterval sourceBlock = Views.offsetInterval(source, sourceMin, sourceSize); + final Img targetBlock = ArrayImgs.unsignedBytes(targetSize); + Downsample.downsample(sourceBlock, targetBlock, factors); - // back up the original block (if present) before overwriting. - final DataBlock originalBlock = reader.readBlock(toDataset, toAttributes, gridPosition); - if (originalBlock != null) { - backupWriter.writeBlock(toDataset, toAttributes, originalBlock); - } - - // delete first so a block that became empty does not leave a stale remnant. - N5Utils.deleteBlock(targetBlock, writer, toDataset, gridPosition); - N5Utils.saveNonEmptyBlock(targetBlock, writer, toDataset, gridPosition, new UnsignedByteType()); - LOG.info("downsampleBlock: updated {} block {}", toDataset, Arrays.toString(gridPosition)); + // back up the original block (if present) before overwriting. + final DataBlock originalBlock = reader.readBlock(toDataset, toAttributes, gridPosition); + if (originalBlock != null) { + backupWriter.writeBlock(toDataset, toAttributes, originalBlock); } + + // delete first so a block that became empty does not leave a stale remnant. + N5Utils.deleteBlock(targetBlock, writer, toDataset, gridPosition); + N5Utils.saveNonEmptyBlock(targetBlock, writer, toDataset, gridPosition, new UnsignedByteType()); + LOG.info("downsampleBlock: updated {} block {}", toDataset, Arrays.toString(gridPosition)); } // ------------------------------------------------------------------------------------------------ From c021df39ec660499c7e57af549bc75d44b1f990a Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sun, 12 Jul 2026 17:03:11 -0400 Subject: [PATCH 10/15] Get sfov positions from render --- .../alignment/multisem/MultiSemUtilities.java | 30 ++ .../multisem/MultiSemUtilitiesTest.java | 11 + .../render/client/TileReorderingClient.java | 19 +- .../spark/multisem/Wafer6061Inpainter.java | 400 +++++++++--------- 4 files changed, 250 insertions(+), 210 deletions(-) diff --git a/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java b/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java index 7bca85591..919461b3d 100644 --- a/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java +++ b/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java @@ -263,6 +263,36 @@ public static boolean isSimpleMFOVName(final String name) { /** Each MFOV has 91 SFOVs or tiles */ public static int NUMBER_OF_TILES_IN_MFOV = 91; + /** + * Maps an SFOV's original spiral acquisition number (the 1-based {@code _s##} value in a multi-SEM tileId, + * assigned by spiraling counterclockwise out from the center of the MFOV) to its 1-based row-major beam index + * (numbering the 91 beams of the MFOV from 1 in the top-left corner, going row by row). This "spiral -> + * row-major" permutation of the 91-beam hexagonal layout is the single source of truth shared by + * {@code TileReorderingClient} (which uses it to order tiles) and the multi-SEM inpainter (which uses it to + * match render tiles to the acquisition xlog's row-major sfov axis). + * + * @param spiralSFOVNumber the 1-based spiral sfov number (1..91). + * @return the corresponding 1-based row-major beam index. + */ + public static int getRowMajorSFOVIndex(final int spiralSFOVNumber) { + return SPIRAL_TO_ROW_MAJOR_SFOV[spiralSFOVNumber - 1]; + } + + // Number the 91 beams of an MFOV from 1 to 91 in row-major order (top-left corner, row by row), then record + // those numbers in the original spiral acquisition order (center-out, counterclockwise). Thus + // SPIRAL_TO_ROW_MAJOR_SFOV[spiralNumber - 1] is that beam's row-major index. + private static final int[] SPIRAL_TO_ROW_MAJOR_SFOV = { + 46, 47, 36, 35, 45, 56, 57, 48, 37, 27, + 26, 25, 34, 44, 55, 65, 66, 67, 58, 49, + 38, 28, 19, 18, 17, 16, 24, 33, 43, 54, + 64, 73, 74, 75, 76, 68, 59, 50, 39, 29, + 20, 12, 11, 10, 9, 8, 15, 23, 32, 42, + 53, 63, 72, 80, 81, 82, 83, 84, 77, 69, + 60, 51, 40, 30, 21, 13, 6, 5, 4, 3, + 2, 1, 7, 14, 22, 31, 41, 52, 62, 71, + 79, 86, 87, 88, 89, 90, 91, 85, 78, 70, 61 + }; + private static final Pattern SIMPLE_MFOV_NAME_PATTERN = Pattern.compile("^m(\\d{4})$"); private static final Logger LOG = LoggerFactory.getLogger(MultiSemUtilities.class); diff --git a/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java b/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java index 4ee758f5a..6889f035a 100644 --- a/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java +++ b/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java @@ -33,4 +33,15 @@ public void testTileIdParsers() { "16", MultiSemUtilities.getSFOVIndexForTileId(manyScanTileId)); } + @Test + public void testGetRowMajorSFOVIndex() { + // spiral -> row-major permutation of the 91-beam hex layout (see TileReorderingClient's original table) + Assert.assertEquals("spiral 1 (center) should be row-major 46", + 46, MultiSemUtilities.getRowMajorSFOVIndex(1)); + Assert.assertEquals("spiral 72 should be row-major 1 (top-left)", + 1, MultiSemUtilities.getRowMajorSFOVIndex(72)); + Assert.assertEquals("spiral 91 should be row-major 61", + 61, MultiSemUtilities.getRowMajorSFOVIndex(91)); + } + } diff --git a/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java b/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java index bb94ba2da..7234bd40b 100644 --- a/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java +++ b/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java @@ -2,6 +2,7 @@ import com.beust.jcommander.Parameter; import com.beust.jcommander.ParametersDelegate; +import org.janelia.alignment.multisem.MultiSemUtilities; import org.janelia.alignment.spec.ResolvedTileSpecCollection; import org.janelia.alignment.spec.TileSpec; import org.janelia.alignment.spec.stack.StackMetaData; @@ -148,22 +149,6 @@ public enum RenderingOrder implements Comparator { }); - // The new order of the tiles in the multi-sem stack: - // Number the tiles in an SFOV from 1 to 91, starting in the top left corner and going - // row by row. Then, record these numbers in the original order of the tiles (i.e., - // starting in the middle of the SFOV and spiraling outwards counterclockwise). - private static final int[] newNumber = { - 46, 47, 36, 35, 45, 56, 57, 48, 37, 27, - 26, 25, 34, 44, 55, 65, 66, 67, 58, 49, - 38, 28, 19, 18, 17, 16, 24, 33, 43, 54, - 64, 73, 74, 75, 76, 68, 59, 50, 39, 29, - 20, 12, 11, 10, 9, 8, 15, 23, 32, 42, - 53, 63, 72, 80, 81, 82, 83, 84, 77, 69, - 60, 51, 40, 30, 21, 13, 6, 5, 4, 3, - 2, 1, 7, 14, 22, 31, 41, 52, 62, 71, - 79, 86, 87, 88, 89, 90, 91, 85, 78, 70, 61 - }; - private static final Pattern TILE_ID_SEPARATOR = Pattern.compile("_"); private final Comparator tileSpecComparator; @@ -196,7 +181,7 @@ private static int reverse(final int order) { } private static int linearIndex(final int sFov) { - return newNumber[sFov - 1]; + return MultiSemUtilities.getRowMajorSFOVIndex(sFov); } private static double[] getUpperEdgeMidpoint(final TileSpec tileSpec) { diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index c2461a69f..a63c7f3e1 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -37,8 +37,12 @@ import org.apache.spark.api.java.JavaRDD; import org.apache.spark.api.java.JavaSparkContext; import org.apache.spark.broadcast.Broadcast; +import org.janelia.alignment.multisem.MultiSemUtilities; +import org.janelia.alignment.spec.Bounds; +import org.janelia.alignment.spec.TileBounds; import org.janelia.alignment.util.Grid; import org.janelia.render.client.ClientRunner; +import org.janelia.render.client.RenderDataClient; import org.janelia.render.client.parameter.CommandLineParameters; import org.janelia.render.client.spark.LogUtilities; import org.janelia.saalfeldlab.n5.DataBlock; @@ -56,10 +60,13 @@ *

* Inputs are (a) an N5 tissue volume (only non-empty blocks stored), (b) an N5 mask marking where image data is * present ({@code mask > 0}) vs. missing ({@code mask == 0}), and (c) the acquisition xlog zarr. The ROI point cloud - * is built in the tissue's s0 voxel frame directly from the xlog: each SFOV's acquisition position ({@code x}, - * {@code y}) is placed exactly as {@code msem_to_render.py} ingests it (rotate by {@code 180 + rotation_slab}, then - * {@code stage - min + halfSFOV}), and carries that SFOV's {@code distance_roi}. No render service and no coordinate - * fitting are needed (x/y/distance_roi share the same {@code (mfov, sfov)} indexing). + * is built in the tissue's s0 voxel frame from the ALIGNED render stack (so montage stitching is accounted for): the + * first layer's tile bounds are fetched from the render web service on the driver, and each tile center (render world + * pixels) is mapped to the voxel frame with {@code voxel = center - translate} (the neuroglancer group-level + * {@code translate}, i.e. the stack bounding-box min). Each tile carries the {@code distance_roi} of its SFOV, read + * from the xlog for the slab whose {@code id_serial} equals {@code serial}: the tile's mfov (0-based) indexes the xlog + * mfov axis directly, and its spiral sfov number (the {@code _s##} field) is mapped to the xlog's row-major sfov axis + * via {@link MultiSemUtilities#getRowMajorSFOVIndex}. *

* The client processes the full-resolution ({@code s0}) blocks of the tissue in parallel. For each block it first * applies the cheap ROI-distance filter (interpolating {@code distance_roi} at the block center via inverse-distance @@ -72,33 +79,20 @@ */ public class Wafer6061Inpainter { - // SFOV image size fallback (pixels), used only when the xlog lacks x_sfov / y_sfov; the xlog values are - // preferred at runtime (see loadPointCloud), so these are wafer 60/61 defaults, not a tuning knob. - private static final int DEFAULT_SFOV_WIDTH = 2000; - private static final int DEFAULT_SFOV_HEIGHT = 1748; - // Inverse-distance-weighting knobs for the ROI-distance interpolation; fixed for wafer 60/61. private static final int IDW_K = 8; private static final double IDW_POWER = 2.0; - // Tissue containers are named w_s_r (e.g. w61_s109_r00); the serial is parsed from that. + // Tissue containers/datasets are named w_s_r (e.g. w61_s109_r00); the serial is parsed from that. private static final Pattern SERIAL_IN_NAME = Pattern.compile("_s(\\d+)"); // xlog field layout (consumed by loadPointCloud). The xlog is a zarr: its arrays are stored C-order, but the // n5-zarr reader REVERSES the axis order, so the imglib2 axis indices below are the reverse of the C-order shape. - // field C-order shape imglib2 axes (what this code sees) meaning - // id_serial (slab,) [slab] serial label per slab position - // x, y (scan, slab, mfov, sfov) [sfov, mfov, slab, scan] SFOV centre, full-res px - // rotation_slab (scan, slab) [slab, scan] per-slab rotation, degrees - // distance_roi (slab, mfov, sfov) [sfov, mfov, slab] distance to ROI, um (scan-independent) - // x_sfov, y_sfov (width,), (height,) [size] SFOV image size, px - // The slab/scan axis indices are HARDCODED (below) rather than discovered by matching sizes: the sizes are not - // guaranteed to be unique. x / y / distance_roi share (mfov, sfov) indexing, so after slicing the slab (and scan) - // axis they all reduce to the same 2-D (sfov, mfov) plane with matching indices. - private static final int XY_AXIS_SLAB = 2; - private static final int XY_AXIS_SCAN = 3; - private static final int ROT_AXIS_SLAB = 0; - private static final int ROT_AXIS_SCAN = 1; + // field C-order shape imglib2 axes (what this code sees) meaning + // id_serial (slab,) [slab] serial label per slab position + // distance_roi (slab, mfov, sfov) [sfov, mfov, slab] distance to ROI, um (scan-independent) + // The slab axis index is HARDCODED (below) rather than discovered by matching sizes (sizes are not guaranteed to + // be unique). After slicing the slab axis, distance_roi reduces to a 2-D (sfov, mfov) plane. private static final int DIST_AXIS_SLAB = 2; private static final Logger LOG = LoggerFactory.getLogger(Wafer6061Inpainter.class); @@ -128,10 +122,33 @@ public static class Parameters extends CommandLineParameters { required = true) public String xlogPath; + @Parameter( + names = "--baseDataUrl", + description = "Base render web service URL for data (e.g. http://host[:port]/render-ws/v1); the aligned " + + "tile positions that define the ROI cloud are read from here.", + required = true) + public String baseDataUrl; + + @Parameter( + names = "--owner", + description = "Render stack owner.") + public String owner = "hess_wafers_60_61"; + + @Parameter( + names = "--project", + description = "Render project holding the aligned stack. Defaults to the name of the directory that " + + "contains the N5 container (e.g. .../w61_serial_100_to_109/w61_s109_r00.n5 -> w61_serial_100_to_109).") + public String project; + + @Parameter( + names = "--stack", + description = "Render stack whose first layer's tile positions define the ROI cloud. Defaults to the dataset name.") + public String stack; + @Parameter( names = "--serial", description = "Serial label (id_serial) of the section stored in this N5, resolved to the reference " + - "arrays' slab position. Defaults to the serial parsed from the container name " + + "arrays' slab position. Defaults to the serial parsed from the dataset name " + "(w_s_r, e.g. w61_s109_r00 -> 109).") public long serial = -1; @@ -173,26 +190,52 @@ public void validate() { throw new IllegalArgumentException("--maxRoiDistance must be positive"); } if (serial < 0) { - serial = inferSerial(n5Path); + serial = inferSerial(dataset); } - } - - /** Parses the serial label from a container name like {@code .../w61_s109_r00} (basename {@code _s}). */ - static long inferSerial(final String n5Path) { - String basename = n5Path; - while (basename.endsWith("/")) { - basename = basename.substring(0, basename.length() - 1); + if (project == null) { + project = inferProject(n5Path); } - final int slash = basename.lastIndexOf('/'); - if (slash >= 0) { - basename = basename.substring(slash + 1); + if (stack == null) { + stack = basename(dataset); } - final Matcher matcher = SERIAL_IN_NAME.matcher(basename); + } + + /** Parses the serial label from a dataset/container name like {@code .../w61_s109_r00} (basename {@code _s}). */ + static long inferSerial(final String path) { + final Matcher matcher = SERIAL_IN_NAME.matcher(basename(path)); if (matcher.find()) { return Long.parseLong(matcher.group(1)); } - throw new IllegalArgumentException("could not infer the serial from n5Path '" + n5Path + - "'; pass --serial explicitly"); + throw new IllegalArgumentException("could not infer the serial from '" + path + "'; pass --serial explicitly"); + } + + /** Infers the render project from the name of the directory that contains the N5 container. */ + static String inferProject(final String n5Path) { + String p = n5Path; + while (p.endsWith("/")) { + p = p.substring(0, p.length() - 1); + } + final int containerSlash = p.lastIndexOf('/'); + if (containerSlash <= 0) { + throw new IllegalArgumentException("could not infer the render project from n5Path '" + n5Path + + "' (no parent directory); pass --project explicitly"); + } + final String project = basename(p.substring(0, containerSlash)); + if (project.isEmpty()) { + throw new IllegalArgumentException("could not infer the render project from n5Path '" + n5Path + + "'; pass --project explicitly"); + } + return project; + } + + /** Last path segment of a path, stripped of any leading group/parent path and trailing slashes. */ + static String basename(final String path) { + String p = path; + while (p.endsWith("/")) { + p = p.substring(0, p.length() - 1); + } + final int slash = p.lastIndexOf('/'); + return slash >= 0 ? p.substring(slash + 1) : p; } } @@ -222,13 +265,8 @@ public void runClient(final String[] args) throws Exception { public void run() throws IOException { - // 1. Load the (small) per-slab 2-D point cloud from the xlog on the driver, placed in the tissue voxel frame. - final PointCloud cloud = loadPointCloud(params.xlogPath, params.serial); - LOG.info("run: loaded {} ROI reference points for serial {}", cloud.size(), params.serial); - - // Read the tissue s0 metadata and discover the existing pyramid levels. - final long[] tissueTranslate; - final long[] maskTranslate; + // Read the tissue s0 metadata, the world->voxel offset, and discover the existing pyramid levels. + final double[] worldToVoxel; final int numDimensions; final List levels = new ArrayList<>(); final Map levelAttributes = new LinkedHashMap<>(); @@ -268,12 +306,10 @@ public void run() throws IOException { Arrays.toString(s0Attributes.getBlockSize())); } - tissueTranslate = readTranslate(n5, params.fullDataset(), numDimensions); - maskTranslate = readTranslate(n5, params.mask, numDimensions); - if (! Arrays.equals(tissueTranslate, maskTranslate)) { - throw new IllegalArgumentException("tissue translate " + Arrays.toString(tissueTranslate) + - " must match mask translate " + Arrays.toString(maskTranslate)); - } + // World->voxel offset for placing render tile centers: the neuroglancer 'translate' (the stack + // bounding-box min in world pixels) is written on the multiscale GROUP, not on s0 (s0 only carries a + // sub-pixel centering transform). May be null here; loadPointCloud falls back to the render stack bounds. + worldToVoxel = readGroupTranslate(n5, params.dataset); levels.add(params.fullDataset()); levelAttributes.put(params.fullDataset(), s0Attributes); @@ -293,9 +329,14 @@ public void run() throws IOException { downsampleFactors = readDownsamplingFactors(n5, levels.get(1), numDimensions); } } - LOG.info("run: tissue translate={}, mask translate={}, pyramid levels={}, downsampleFactors={}", - Arrays.toString(tissueTranslate), Arrays.toString(maskTranslate), levels, - Arrays.toString(downsampleFactors)); + LOG.info("run: world->voxel translate={}, pyramid levels={}, downsampleFactors={}", + Arrays.toString(worldToVoxel), levels, Arrays.toString(downsampleFactors)); + + // Load the (small) per-slab ROI point cloud on the driver: positions from the aligned render stack, distances + // from the xlog. A single render request keeps the server load light; the cloud is then broadcast to executors. + final PointCloud cloud = loadPointCloud(params, worldToVoxel); + LOG.info("run: loaded {} ROI reference points for serial {} (render {}/{}/{})", + cloud.size(), params.serial, params.owner, params.project, params.stack); // Create the backup container and mirror all datasets that might receive backups. final String backupPath = params.getBackupPath(); @@ -370,142 +411,122 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, } // ------------------------------------------------------------------------------------------------ - // Step 1: load the per-slab 2-D point cloud from the xlog zarr. + // Step 1: load the per-slab ROI point cloud (positions from the aligned render stack, distances from the xlog). // ------------------------------------------------------------------------------------------------ /** - * Builds the per-slab ROI point cloud in the tissue s0 voxel frame, entirely from the xlog. For the slab - * whose {@code id_serial} equals {@code serial}, each SFOV's acquisition position ({@code x}, {@code y}, in - * full-resolution pixels) is placed into the render/ingestion frame exactly as {@code msem_to_render.py} does: - * 

-	 *   center = EuclideanTransform(rotation = radians(180 + rotation_slab)) . (x, y)
-	 *   voxel  = center - min(center) + (sfovWidth/2, sfovHeight/2)
-	 *
- * i.e. the ingestion {@code stage - min + margin} placement (the constant {@code margin} and the export - * {@code translate} cancel out into the voxel frame). Each voxel carries the scan-independent {@code distance_roi} - * for that SFOV. Since {@code x}, {@code y} and {@code distance_roi} are all indexed by {@code (mfov, sfov)} in the - * xlog, the correspondence is exact and no fitting is needed. Alignment (montage stitching) only perturbs these - * positions slightly, well within {@code maxRoiDistance}, so the unaligned ingestion placement is used directly. + * Builds the per-slab ROI point cloud in the tissue s0 voxel frame. Positions come from the ALIGNED render + * stack (so montage stitching is accounted for): the first layer's tile bounds are fetched from the render web + * service on the driver, and each tile center (render world pixels) is mapped to the voxel frame with + * {@code voxel = center - worldToVoxel} ({@code worldToVoxel} is the neuroglancer group-level {@code translate}, + * i.e. the stack bounding-box min; if it is null we fall back to the render stack bounds). Each tile carries the + * {@code distance_roi} of its SFOV, read from the xlog for the slab whose {@code id_serial} equals + * {@code params.serial}: the tile's mfov (0-based, parsed from the tileId) indexes the xlog mfov axis directly, and + * its spiral sfov number (the {@code _s##} field) is mapped to the xlog's row-major sfov axis via + * {@link MultiSemUtilities#getRowMajorSFOVIndex}. Tiles whose SFOV has no (NaN) {@code distance_roi}, or an + * out-of-range mfov/sfov, are dropped. */ - static PointCloud loadPointCloud(final String xlogPath, - final long serial) { - try (final N5Reader xlog = new N5Factory().openReader(xlogPath)) { + static PointCloud loadPointCloud(final Parameters params, + final double[] worldToVoxel) + throws IOException { + + // (a) xlog: the slab's 2-D distance_roi grid, materialized into a plain array so it outlives the xlog reader. + final double[][] distBySfovMfov; // [rowMajorSfov0Based][mfov0Based] + try (final N5Reader xlog = new N5Factory().openReader(params.xlogPath)) { final double[] idSerial = read1d(xlog, "id_serial"); - final int slabPosition = findSlabPosition(idSerial, serial); + final int slabPosition = findSlabPosition(idSerial, params.serial); if (slabPosition < 0) { final long[] available = new long[idSerial.length]; for (int i = 0; i < idSerial.length; i++) { available[i] = Math.round(idSerial[i]); } - throw new IllegalArgumentException("serial " + serial + " not found in id_serial; available serials are " + + throw new IllegalArgumentException("serial " + params.serial + " not found in id_serial; available serials are " + Arrays.toString(available)); } - final long slabCount = idSerial.length; // 413 for wafer 61 - // SFOV image size defines the half-SFOV placement offset; prefer the xlog (x_sfov / y_sfov), fall back to constants. - final int sfW = sfovSize(xlog, "x_sfov", DEFAULT_SFOV_WIDTH); - final int sfH = sfovSize(xlog, "y_sfov", DEFAULT_SFOV_HEIGHT); - - // Open the reference arrays (axis layout documented and hardcoded at the top of the class) and verify each - // one carries the slab count on its expected axis before we slice by it. - final RandomAccessibleInterval xAll = openDoubles(xlog, "x"); - final RandomAccessibleInterval yAll = openDoubles(xlog, "y"); - final RandomAccessibleInterval rotAll = openDoubles(xlog, "rotation_slab"); final RandomAccessibleInterval distAll = openDoubles(xlog, "distance_roi"); - requireSlabAxis(xAll, XY_AXIS_SLAB, slabCount, "x"); - requireSlabAxis(yAll, XY_AXIS_SLAB, slabCount, "y"); - requireSlabAxis(rotAll, ROT_AXIS_SLAB, slabCount, "rotation_slab"); - requireSlabAxis(distAll, DIST_AXIS_SLAB, slabCount, "distance_roi"); - final long nScans = rotAll.dimension(ROT_AXIS_SCAN); + requireSlabAxis(distAll, DIST_AXIS_SLAB, idSerial.length, "distance_roi"); // distance_roi is scan-independent: slice the slab axis -> 2-D (sfov, mfov). final RandomAccessibleInterval distSlab = Views.hyperSlice(distAll, DIST_AXIS_SLAB, slabPosition); - - // Choose the scan whose x/y define the cloud: the first with finite x and rotation_slab for this slab. - // The positions are nearly scan-independent, so no scan override is needed. - int scan = -1; - for (int s = 0; s < nScans; s++) { - if (! Double.isFinite(rotationAt(rotAll, slabPosition, s))) { - continue; - } - if (hasFinite(sliceScanAndSlab(xAll, s, slabPosition))) { - scan = s; - break; - } - } - if (scan < 0) { - throw new IllegalArgumentException("no scan with finite x and rotation_slab found for serial " + - serial + " (slab position " + slabPosition + ")"); - } - - final double rotationSlab = rotationAt(rotAll, slabPosition, scan); - final double theta = Math.toRadians(180.0 + rotationSlab); - final double cos = Math.cos(theta); - final double sin = Math.sin(theta); - - final RandomAccessibleInterval xSlab = sliceScanAndSlab(xAll, scan, slabPosition); - final RandomAccessibleInterval ySlab = sliceScanAndSlab(yAll, scan, slabPosition); - - // Place each SFOV center (rotation only) and keep the finite ones with their distance_roi. - final List cxs = new ArrayList<>(); - final List cys = new ArrayList<>(); - final List dists = new ArrayList<>(); - final RandomAccess xra = xSlab.randomAccess(); - final RandomAccess yra = ySlab.randomAccess(); + final int nSfov = (int) distSlab.dimension(0); + final int nMfov = (int) distSlab.dimension(1); + distBySfovMfov = new double[nSfov][nMfov]; final RandomAccess dra = distSlab.randomAccess(); final long[] pos = new long[2]; - double minX = Double.POSITIVE_INFINITY; - double minY = Double.POSITIVE_INFINITY; - for (long i0 = 0; i0 < xSlab.dimension(0); i0++) { - for (long i1 = 0; i1 < xSlab.dimension(1); i1++) { - pos[0] = i0; - pos[1] = i1; - final double x = xra.setPositionAndGet(pos).get(); - final double y = yra.setPositionAndGet(pos).get(); - final double d = dra.setPositionAndGet(pos).get(); - if (Double.isFinite(x) && Double.isFinite(y) && Double.isFinite(d)) { - final double cx = cos * x - sin * y; - final double cy = sin * x + cos * y; - cxs.add(cx); - cys.add(cy); - dists.add(d); - minX = Math.min(minX, cx); - minY = Math.min(minY, cy); - } + for (int s = 0; s < nSfov; s++) { + for (int m = 0; m < nMfov; m++) { + pos[0] = s; + pos[1] = m; + distBySfovMfov[s][m] = dra.setPositionAndGet(pos).get(); } } - if (cxs.isEmpty()) { - throw new IllegalArgumentException("no finite ROI points found for serial " + serial + - " (slab position " + slabPosition + ", scan " + scan + ")"); - } - - // Shift into the voxel frame: min(center) -> half-SFOV (so the min tile's top-left is at the origin). - final double[] xs = new double[cxs.size()]; - final double[] ys = new double[cys.size()]; - final double[] ds = new double[dists.size()]; - for (int i = 0; i < xs.length; i++) { - xs[i] = cxs.get(i) - minX + sfW / 2.0; - ys[i] = cys.get(i) - minY + sfH / 2.0; - ds[i] = dists.get(i); - } - LOG.info("loadPointCloud: serial {} -> slab position {}, scan {}, rotation_slab {} deg, sfov {}x{}, {} points", - serial, slabPosition, scan, rotationSlab, sfW, sfH, xs.length); - return new PointCloud(xs, ys, ds); + LOG.info("loadPointCloud: serial {} -> slab position {}, distance_roi grid is {} sfov x {} mfov", + params.serial, slabPosition, nSfov, nMfov); + } + final int nSfov = distBySfovMfov.length; + final int nMfov = nSfov > 0 ? distBySfovMfov[0].length : 0; + + // (b) render: the first layer's aligned tile centers, fetched once on the driver. + final RenderDataClient renderClient = new RenderDataClient(params.baseDataUrl, params.owner, params.project); + final List zValues = renderClient.getStackZValues(params.stack); + if (zValues.isEmpty()) { + throw new IllegalArgumentException("render stack " + params.owner + "/" + params.project + "/" + + params.stack + " has no z layers"); + } + final double firstZ = zValues.get(0); + final List tiles = renderClient.getTileBounds(params.stack, firstZ); + + final double offsetX; + final double offsetY; + if (worldToVoxel != null) { + offsetX = worldToVoxel[0]; + offsetY = worldToVoxel[1]; + } else { + final Bounds stackBounds = renderClient.getStackMetaData(params.stack).getStackBounds(); + if (stackBounds == null || stackBounds.getMinX() == null || stackBounds.getMinY() == null) { + throw new IllegalArgumentException("N5 group " + params.dataset + " has no neuroglancer 'translate' and " + + "render stack " + params.stack + " has no bounds; cannot map world coordinates to voxels"); + } + offsetX = stackBounds.getMinX(); + offsetY = stackBounds.getMinY(); + LOG.warn("loadPointCloud: N5 group {} has no 'translate'; using render stack bounds min ({}, {}) as the world->voxel offset", + params.dataset, offsetX, offsetY); } - } - /** SFOV image size from the xlog {@code x_sfov}/{@code y_sfov} length, or {@code fallback} if that dataset is absent. */ - private static int sfovSize(final N5Reader xlog, final String dataset, final int fallback) { - try { - final DatasetAttributes attributes = xlog.getDatasetAttributes(dataset); - if (attributes != null && attributes.getNumDimensions() >= 1) { - return (int) attributes.getDimensions()[0]; - } - } catch (final Exception e) { - LOG.warn("sfovSize: could not read {} ({}), using fallback {}", dataset, e.getMessage(), fallback); + // Attach each tile's distance_roi (by mfov + spiral->row-major sfov) and place its center in the voxel frame. + final List xs = new ArrayList<>(); + final List ys = new ArrayList<>(); + final List ds = new ArrayList<>(); + int noDistance = 0; + int outOfRange = 0; + for (final TileBounds tile : tiles) { + final String tileId = tile.getTileId(); + final int mfov = Integer.parseInt(MultiSemUtilities.getSimpleMfovForTileId(tileId).substring(1)); // m0013 -> 13 + final int spiralSfov = Integer.parseInt(MultiSemUtilities.getSFOVIndexForTileId(tileId)); // _s## spiral, 1-based + final int sfov = MultiSemUtilities.getRowMajorSFOVIndex(spiralSfov) - 1; // xlog row-major, 0-based + if (mfov < 0 || mfov >= nMfov || sfov < 0 || sfov >= nSfov) { + outOfRange++; + continue; + } + final double d = distBySfovMfov[sfov][mfov]; + if (! Double.isFinite(d)) { + noDistance++; + continue; + } + xs.add(tile.getCenterX() - offsetX); + ys.add(tile.getCenterY() - offsetY); + ds.add(d); } - return fallback; + if (xs.isEmpty()) { + throw new IllegalArgumentException("no render tiles in stack " + params.stack + " z " + firstZ + + " could be matched to a finite xlog distance_roi (fetched " + tiles.size() + + " tiles; " + noDistance + " had NaN distance, " + outOfRange + " had out-of-range mfov/sfov)"); + } + LOG.info("loadPointCloud: matched {} of {} render tiles (z {}) to distance_roi; dropped {} NaN-distance, {} out-of-range; world->voxel offset ({}, {})", + xs.size(), tiles.size(), firstZ, noDistance, outOfRange, offsetX, offsetY); + return new PointCloud(toArray(xs), toArray(ys), toArray(ds)); } /** Fails fast if a hardcoded slab axis does not carry the expected slab count (guards the layout at the top). */ @@ -519,32 +540,13 @@ private static void requireSlabAxis(final RandomAccessibleInterval img, final } } - /** Value of the 2-D {@code rotation_slab} (axes {@code [slab, scan]}) at the given slab position and scan. */ - private static double rotationAt(final RandomAccessibleInterval rotAll, - final long slabPosition, - final long scan) { - final long[] p = new long[2]; - p[ROT_AXIS_SLAB] = slabPosition; - p[ROT_AXIS_SCAN] = scan; - return rotAll.randomAccess().setPositionAndGet(p).get(); - } - - /** Slices the 4-D {@code x}/{@code y} array (axes {@code [sfov, mfov, slab, scan]}) to the 2-D (sfov, mfov) plane. */ - private static RandomAccessibleInterval sliceScanAndSlab(final RandomAccessibleInterval xy, - final long scan, - final long slabPosition) { - // slice the higher axis (scan) first so the lower slab axis index stays valid. - return Views.hyperSlice(Views.hyperSlice(xy, XY_AXIS_SCAN, scan), XY_AXIS_SLAB, slabPosition); - } - - /** True if the interval has at least one finite value. */ - private static boolean hasFinite(final RandomAccessibleInterval img) { - for (final DoubleType t : Views.iterable(img)) { - if (Double.isFinite(t.get())) { - return true; - } + /** Copies a list of doubles into a primitive array. */ + private static double[] toArray(final List list) { + final double[] array = new double[list.size()]; + for (int i = 0; i < array.length; i++) { + array[i] = list.get(i); } - return false; + return array; } /** Reads a 1-D double dataset (blosc-safe: uses the block-not-found overload to avoid the getAttribute NPE). */ @@ -878,9 +880,21 @@ static long[] affectedBlock(final long[] previousGridPosition, final int[] facto return g; } - private static long[] readTranslate(final N5Reader n5, final String dataset, final int numDimensions) { - final long[] translate = n5.getAttribute(dataset, "translate", long[].class); - return translate != null ? translate : new long[numDimensions]; + /** + * Reads the neuroglancer {@code translate} (the stack bounding-box min in world pixels) from the multiscale group, + * or returns null if it is absent. The render N5 export writes it on the group ({@code }), not on + * {@code s0} (s0 only carries a sub-pixel centering {@code transform}); callers fall back to the render stack bounds. + */ + private static double[] readGroupTranslate(final N5Reader n5, final String group) { + try { + final double[] translate = n5.getAttribute(group, "translate", double[].class); + if (translate != null && translate.length >= 2) { + return translate; + } + } catch (final Exception e) { + LOG.warn("readGroupTranslate: could not read 'translate' from group {} ({})", group, e.getMessage()); + } + return null; } /** Opens an N5 reader for a tissue/backup container (explicit N5 format; local path or gs://). */ From 49cb4c5c5124b6fcbcd0df006dc4d67a3b4895b9 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Sun, 12 Jul 2026 17:21:06 -0400 Subject: [PATCH 11/15] Read render info from N5 attributes --- .../spark/multisem/Wafer6061Inpainter.java | 170 ++++++++++-------- 1 file changed, 96 insertions(+), 74 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index a63c7f3e1..a7d298900 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -122,29 +122,6 @@ public static class Parameters extends CommandLineParameters { required = true) public String xlogPath; - @Parameter( - names = "--baseDataUrl", - description = "Base render web service URL for data (e.g. http://host[:port]/render-ws/v1); the aligned " + - "tile positions that define the ROI cloud are read from here.", - required = true) - public String baseDataUrl; - - @Parameter( - names = "--owner", - description = "Render stack owner.") - public String owner = "hess_wafers_60_61"; - - @Parameter( - names = "--project", - description = "Render project holding the aligned stack. Defaults to the name of the directory that " + - "contains the N5 container (e.g. .../w61_serial_100_to_109/w61_s109_r00.n5 -> w61_serial_100_to_109).") - public String project; - - @Parameter( - names = "--stack", - description = "Render stack whose first layer's tile positions define the ROI cloud. Defaults to the dataset name.") - public String stack; - @Parameter( names = "--serial", description = "Serial label (id_serial) of the section stored in this N5, resolved to the reference " + @@ -192,12 +169,6 @@ public void validate() { if (serial < 0) { serial = inferSerial(dataset); } - if (project == null) { - project = inferProject(n5Path); - } - if (stack == null) { - stack = basename(dataset); - } } /** Parses the serial label from a dataset/container name like {@code .../w61_s109_r00} (basename {@code _s}). */ @@ -209,25 +180,6 @@ static long inferSerial(final String path) { throw new IllegalArgumentException("could not infer the serial from '" + path + "'; pass --serial explicitly"); } - /** Infers the render project from the name of the directory that contains the N5 container. */ - static String inferProject(final String n5Path) { - String p = n5Path; - while (p.endsWith("/")) { - p = p.substring(0, p.length() - 1); - } - final int containerSlash = p.lastIndexOf('/'); - if (containerSlash <= 0) { - throw new IllegalArgumentException("could not infer the render project from n5Path '" + n5Path + - "' (no parent directory); pass --project explicitly"); - } - final String project = basename(p.substring(0, containerSlash)); - if (project.isEmpty()) { - throw new IllegalArgumentException("could not infer the render project from n5Path '" + n5Path + - "'; pass --project explicitly"); - } - return project; - } - /** Last path segment of a path, stripped of any leading group/parent path and trailing slashes. */ static String basename(final String path) { String p = path; @@ -265,8 +217,9 @@ public void runClient(final String[] args) throws Exception { public void run() throws IOException { - // Read the tissue s0 metadata, the world->voxel offset, and discover the existing pyramid levels. + // Read the tissue s0 metadata, the render target, the world->voxel offset, and discover the pyramid levels. final double[] worldToVoxel; + final RenderTarget renderTarget; final int numDimensions; final List levels = new ArrayList<>(); final Map levelAttributes = new LinkedHashMap<>(); @@ -306,6 +259,11 @@ public void run() throws IOException { Arrays.toString(s0Attributes.getBlockSize())); } + // Render service parameters (baseDataUrl / owner / project / stack) come straight from the group's + // "renderExport" metadata written by render's N5 export — the same attributes.json that holds the pyramid + // scales and translate — so they never have to be passed on the command line. + renderTarget = readRenderTarget(n5, params.dataset); + // World->voxel offset for placing render tile centers: the neuroglancer 'translate' (the stack // bounding-box min in world pixels) is written on the multiscale GROUP, not on s0 (s0 only carries a // sub-pixel centering transform). May be null here; loadPointCloud falls back to the render stack bounds. @@ -322,11 +280,12 @@ public void run() throws IOException { levelAttributes.put(levelDataset, n5.getDatasetAttributes(levelDataset)); } - // The relative per-level downsampling factor is read from the pyramid itself rather than passed in: s1's - // "downsamplingFactors" attribute is the factor relative to s0, and these pyramids are built with a - // constant factor at every step (see DownsampleHelper / N5DownsamplerSpark), so it applies to all levels. + // The relative per-level downsampling factor is read from the pyramid itself rather than passed in: the + // group's neuroglancer "scales" attribute lists the cumulative factor per level relative to s0, and these + // pyramids are built with a constant factor at every step (see DownsampleHelper / N5DownsamplerSpark), so + // scales[1] applies to all levels. if (levels.size() > 1) { - downsampleFactors = readDownsamplingFactors(n5, levels.get(1), numDimensions); + downsampleFactors = readDownsamplingFactors(n5, params.dataset, numDimensions); } } LOG.info("run: world->voxel translate={}, pyramid levels={}, downsampleFactors={}", @@ -334,9 +293,12 @@ public void run() throws IOException { // Load the (small) per-slab ROI point cloud on the driver: positions from the aligned render stack, distances // from the xlog. A single render request keeps the server load light; the cloud is then broadcast to executors. - final PointCloud cloud = loadPointCloud(params, worldToVoxel); - LOG.info("run: loaded {} ROI reference points for serial {} (render {}/{}/{})", - cloud.size(), params.serial, params.owner, params.project, params.stack); + final RenderDataClient renderClient = + new RenderDataClient(renderTarget.baseDataUrl, renderTarget.owner, renderTarget.project); + final PointCloud cloud = loadPointCloud(params, renderClient, renderTarget.stack, worldToVoxel); + LOG.info("run: loaded {} ROI reference points for serial {} (render {} {}/{}/{})", + cloud.size(), params.serial, renderTarget.baseDataUrl, renderTarget.owner, renderTarget.project, + renderTarget.stack); // Create the backup container and mirror all datasets that might receive backups. final String backupPath = params.getBackupPath(); @@ -427,6 +389,8 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, * out-of-range mfov/sfov, are dropped. */ static PointCloud loadPointCloud(final Parameters params, + final RenderDataClient renderClient, + final String stack, final double[] worldToVoxel) throws IOException { @@ -469,14 +433,12 @@ static PointCloud loadPointCloud(final Parameters params, final int nMfov = nSfov > 0 ? distBySfovMfov[0].length : 0; // (b) render: the first layer's aligned tile centers, fetched once on the driver. - final RenderDataClient renderClient = new RenderDataClient(params.baseDataUrl, params.owner, params.project); - final List zValues = renderClient.getStackZValues(params.stack); + final List zValues = renderClient.getStackZValues(stack); if (zValues.isEmpty()) { - throw new IllegalArgumentException("render stack " + params.owner + "/" + params.project + "/" + - params.stack + " has no z layers"); + throw new IllegalArgumentException("render stack " + stack + " has no z layers"); } final double firstZ = zValues.get(0); - final List tiles = renderClient.getTileBounds(params.stack, firstZ); + final List tiles = renderClient.getTileBounds(stack, firstZ); final double offsetX; final double offsetY; @@ -484,10 +446,10 @@ static PointCloud loadPointCloud(final Parameters params, offsetX = worldToVoxel[0]; offsetY = worldToVoxel[1]; } else { - final Bounds stackBounds = renderClient.getStackMetaData(params.stack).getStackBounds(); + final Bounds stackBounds = renderClient.getStackMetaData(stack).getStackBounds(); if (stackBounds == null || stackBounds.getMinX() == null || stackBounds.getMinY() == null) { throw new IllegalArgumentException("N5 group " + params.dataset + " has no neuroglancer 'translate' and " + - "render stack " + params.stack + " has no bounds; cannot map world coordinates to voxels"); + "render stack " + stack + " has no bounds; cannot map world coordinates to voxels"); } offsetX = stackBounds.getMinX(); offsetY = stackBounds.getMinY(); @@ -505,8 +467,13 @@ static PointCloud loadPointCloud(final Parameters params, final String tileId = tile.getTileId(); final int mfov = Integer.parseInt(MultiSemUtilities.getSimpleMfovForTileId(tileId).substring(1)); // m0013 -> 13 final int spiralSfov = Integer.parseInt(MultiSemUtilities.getSFOVIndexForTileId(tileId)); // _s## spiral, 1-based + // guard the spiral number before the permutation lookup (91-element table) and the mfov before the grid. + if (mfov < 0 || mfov >= nMfov || spiralSfov < 1 || spiralSfov > MultiSemUtilities.NUMBER_OF_TILES_IN_MFOV) { + outOfRange++; + continue; + } final int sfov = MultiSemUtilities.getRowMajorSFOVIndex(spiralSfov) - 1; // xlog row-major, 0-based - if (mfov < 0 || mfov >= nMfov || sfov < 0 || sfov >= nSfov) { + if (sfov < 0 || sfov >= nSfov) { outOfRange++; continue; } @@ -520,7 +487,7 @@ static PointCloud loadPointCloud(final Parameters params, ds.add(d); } if (xs.isEmpty()) { - throw new IllegalArgumentException("no render tiles in stack " + params.stack + " z " + firstZ + + throw new IllegalArgumentException("no render tiles in stack " + stack + " z " + firstZ + " could be matched to a finite xlog distance_roi (fetched " + tiles.size() + " tiles; " + noDistance + " had NaN distance, " + outOfRange + " had out-of-range mfov/sfov)"); } @@ -908,22 +875,77 @@ private static N5Writer openN5Writer(final String path) { } /** - * Reads the {@code downsamplingFactors} attribute of a pyramid level (the factor relative to s0, as written by - * {@code N5DownsamplerSpark} / render's export). Fails fast when it is absent so the pyramid factor is never guessed. + * Reads the per-step pyramid downsampling factor from the multiscale group's neuroglancer {@code scales} attribute + * (as written by render's N5 export). {@code scales} is the cumulative factor per level relative to s0, e.g. + * {@code [[1,1,1],[2,2,1],[4,4,1],...]}, so {@code scales[1]} is the factor of s1 relative to s0 — and because these + * pyramids use a constant step at every level, it is the per-step factor for all levels. Fails fast when it is + * absent or has fewer than two levels, so the pyramid factor is never guessed. */ - private static int[] readDownsamplingFactors(final N5Reader n5, final String dataset, final int numDimensions) { - final int[] factors = n5.getAttribute(dataset, "downsamplingFactors", int[].class); - if (factors == null) { - throw new IllegalArgumentException("dataset " + dataset + " has no 'downsamplingFactors' attribute to derive " + - "the pyramid downsampling factor from"); + private static int[] readDownsamplingFactors(final N5Reader n5, final String group, final int numDimensions) { + final int[][] scales = n5.getAttribute(group, "scales", int[][].class); + if (scales == null || scales.length < 2) { + throw new IllegalArgumentException("group " + group + " has no 'scales' attribute with at least two levels " + + "to derive the pyramid downsampling factor from"); } + final int[] factors = scales[1]; // s1 relative to s0 == the per-step factor (constant-step pyramids) if (factors.length != numDimensions) { - throw new IllegalArgumentException("dataset " + dataset + " has downsamplingFactors " + - Arrays.toString(factors) + " but the volume is " + numDimensions + "-dimensional"); + throw new IllegalArgumentException("group " + group + " has scales[1] " + Arrays.toString(factors) + + " but the volume is " + numDimensions + "-dimensional"); } return factors; } + /** + * Reads the render service coordinates (baseDataUrl / owner / project / stack) straight from the multiscale group's + * {@code renderExport} metadata, written by render's N5 export (the same attributes.json that holds the pyramid + * {@code scales} and {@code translate}). Fails fast when it is absent so the render target is never guessed. + */ + private static RenderTarget readRenderTarget(final N5Reader n5, final String group) { + final RenderExport export = n5.getAttribute(group, "renderExport", RenderExport.class); + if (export == null || export.runParameters == null || export.runParameters.renderWeb == null || + export.runParameters.renderWeb.baseDataUrl == null || export.runParameters.stack == null) { + throw new IllegalArgumentException( + "N5 group " + group + " has no usable 'renderExport' metadata (need runParameters.renderWeb." + + "baseDataUrl/owner/project and runParameters.stack); this client reads the render service parameters from there"); + } + final RenderExport.RenderWeb web = export.runParameters.renderWeb; + return new RenderTarget(web.baseDataUrl, web.owner, web.project, export.runParameters.stack); + } + + /** Render service coordinates resolved from the group's {@code renderExport} metadata. */ + static class RenderTarget { + final String baseDataUrl; + final String owner; + final String project; + final String stack; + + RenderTarget(final String baseDataUrl, final String owner, final String project, final String stack) { + this.baseDataUrl = baseDataUrl; + this.owner = owner; + this.project = project; + this.stack = stack; + } + } + + /** + * Minimal GSON view of the group's {@code renderExport} attribute (only the render coordinates this client needs; + * all other fields written by the export are ignored). + */ + private static class RenderExport { + RunParameters runParameters; + + private static class RunParameters { + RenderWeb renderWeb; + String stack; + } + + private static class RenderWeb { + String baseDataUrl; + String owner; + String project; + } + } + /** Result of inpainting one block: the (block-sized) inpainted image and whether any pixel changed. */ static class InpaintResult { final Img inpainted; From 7797508a7e8dead7e7a85ea8284f5bc9248ef53d Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Mon, 13 Jul 2026 11:49:18 -0400 Subject: [PATCH 12/15] Use correct sfov number --- .../alignment/multisem/MultiSemUtilities.java | 30 ------------------- .../multisem/MultiSemUtilitiesTest.java | 11 ------- .../render/client/TileReorderingClient.java | 19 ++++++++++-- .../spark/multisem/Wafer6061Inpainter.java | 21 +++++-------- 4 files changed, 24 insertions(+), 57 deletions(-) diff --git a/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java b/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java index 919461b3d..7bca85591 100644 --- a/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java +++ b/render-app/src/main/java/org/janelia/alignment/multisem/MultiSemUtilities.java @@ -263,36 +263,6 @@ public static boolean isSimpleMFOVName(final String name) { /** Each MFOV has 91 SFOVs or tiles */ public static int NUMBER_OF_TILES_IN_MFOV = 91; - /** - * Maps an SFOV's original spiral acquisition number (the 1-based {@code _s##} value in a multi-SEM tileId, - * assigned by spiraling counterclockwise out from the center of the MFOV) to its 1-based row-major beam index - * (numbering the 91 beams of the MFOV from 1 in the top-left corner, going row by row). This "spiral -> - * row-major" permutation of the 91-beam hexagonal layout is the single source of truth shared by - * {@code TileReorderingClient} (which uses it to order tiles) and the multi-SEM inpainter (which uses it to - * match render tiles to the acquisition xlog's row-major sfov axis). - * - * @param spiralSFOVNumber the 1-based spiral sfov number (1..91). - * @return the corresponding 1-based row-major beam index. - */ - public static int getRowMajorSFOVIndex(final int spiralSFOVNumber) { - return SPIRAL_TO_ROW_MAJOR_SFOV[spiralSFOVNumber - 1]; - } - - // Number the 91 beams of an MFOV from 1 to 91 in row-major order (top-left corner, row by row), then record - // those numbers in the original spiral acquisition order (center-out, counterclockwise). Thus - // SPIRAL_TO_ROW_MAJOR_SFOV[spiralNumber - 1] is that beam's row-major index. - private static final int[] SPIRAL_TO_ROW_MAJOR_SFOV = { - 46, 47, 36, 35, 45, 56, 57, 48, 37, 27, - 26, 25, 34, 44, 55, 65, 66, 67, 58, 49, - 38, 28, 19, 18, 17, 16, 24, 33, 43, 54, - 64, 73, 74, 75, 76, 68, 59, 50, 39, 29, - 20, 12, 11, 10, 9, 8, 15, 23, 32, 42, - 53, 63, 72, 80, 81, 82, 83, 84, 77, 69, - 60, 51, 40, 30, 21, 13, 6, 5, 4, 3, - 2, 1, 7, 14, 22, 31, 41, 52, 62, 71, - 79, 86, 87, 88, 89, 90, 91, 85, 78, 70, 61 - }; - private static final Pattern SIMPLE_MFOV_NAME_PATTERN = Pattern.compile("^m(\\d{4})$"); private static final Logger LOG = LoggerFactory.getLogger(MultiSemUtilities.class); diff --git a/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java b/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java index 6889f035a..4ee758f5a 100644 --- a/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java +++ b/render-app/src/test/java/org/janelia/alignment/multisem/MultiSemUtilitiesTest.java @@ -33,15 +33,4 @@ public void testTileIdParsers() { "16", MultiSemUtilities.getSFOVIndexForTileId(manyScanTileId)); } - @Test - public void testGetRowMajorSFOVIndex() { - // spiral -> row-major permutation of the 91-beam hex layout (see TileReorderingClient's original table) - Assert.assertEquals("spiral 1 (center) should be row-major 46", - 46, MultiSemUtilities.getRowMajorSFOVIndex(1)); - Assert.assertEquals("spiral 72 should be row-major 1 (top-left)", - 1, MultiSemUtilities.getRowMajorSFOVIndex(72)); - Assert.assertEquals("spiral 91 should be row-major 61", - 61, MultiSemUtilities.getRowMajorSFOVIndex(91)); - } - } diff --git a/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java b/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java index 7234bd40b..bb94ba2da 100644 --- a/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java +++ b/render-ws-java-client/src/main/java/org/janelia/render/client/TileReorderingClient.java @@ -2,7 +2,6 @@ import com.beust.jcommander.Parameter; import com.beust.jcommander.ParametersDelegate; -import org.janelia.alignment.multisem.MultiSemUtilities; import org.janelia.alignment.spec.ResolvedTileSpecCollection; import org.janelia.alignment.spec.TileSpec; import org.janelia.alignment.spec.stack.StackMetaData; @@ -149,6 +148,22 @@ public enum RenderingOrder implements Comparator { }); + // The new order of the tiles in the multi-sem stack: + // Number the tiles in an SFOV from 1 to 91, starting in the top left corner and going + // row by row. Then, record these numbers in the original order of the tiles (i.e., + // starting in the middle of the SFOV and spiraling outwards counterclockwise). + private static final int[] newNumber = { + 46, 47, 36, 35, 45, 56, 57, 48, 37, 27, + 26, 25, 34, 44, 55, 65, 66, 67, 58, 49, + 38, 28, 19, 18, 17, 16, 24, 33, 43, 54, + 64, 73, 74, 75, 76, 68, 59, 50, 39, 29, + 20, 12, 11, 10, 9, 8, 15, 23, 32, 42, + 53, 63, 72, 80, 81, 82, 83, 84, 77, 69, + 60, 51, 40, 30, 21, 13, 6, 5, 4, 3, + 2, 1, 7, 14, 22, 31, 41, 52, 62, 71, + 79, 86, 87, 88, 89, 90, 91, 85, 78, 70, 61 + }; + private static final Pattern TILE_ID_SEPARATOR = Pattern.compile("_"); private final Comparator tileSpecComparator; @@ -181,7 +196,7 @@ private static int reverse(final int order) { } private static int linearIndex(final int sFov) { - return MultiSemUtilities.getRowMajorSFOVIndex(sFov); + return newNumber[sFov - 1]; } private static double[] getUpperEdgeMidpoint(final TileSpec tileSpec) { diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index a7d298900..10721494a 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -65,8 +65,8 @@ * pixels) is mapped to the voxel frame with {@code voxel = center - translate} (the neuroglancer group-level * {@code translate}, i.e. the stack bounding-box min). Each tile carries the {@code distance_roi} of its SFOV, read * from the xlog for the slab whose {@code id_serial} equals {@code serial}: the tile's mfov (0-based) indexes the xlog - * mfov axis directly, and its spiral sfov number (the {@code _s##} field) is mapped to the xlog's row-major sfov axis - * via {@link MultiSemUtilities#getRowMajorSFOVIndex}. + * mfov axis directly, and its sfov number (the 1-based {@code _s##} field) maps directly to the 0-based xlog sfov axis + * as {@code _s## - 1}. *

* The client processes the full-resolution ({@code s0}) blocks of the tissue in parallel. For each block it first * applies the cheap ROI-distance filter (interpolating {@code distance_roi} at the block center via inverse-distance @@ -384,9 +384,8 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, * i.e. the stack bounding-box min; if it is null we fall back to the render stack bounds). Each tile carries the * {@code distance_roi} of its SFOV, read from the xlog for the slab whose {@code id_serial} equals * {@code params.serial}: the tile's mfov (0-based, parsed from the tileId) indexes the xlog mfov axis directly, and - * its spiral sfov number (the {@code _s##} field) is mapped to the xlog's row-major sfov axis via - * {@link MultiSemUtilities#getRowMajorSFOVIndex}. Tiles whose SFOV has no (NaN) {@code distance_roi}, or an - * out-of-range mfov/sfov, are dropped. + * its sfov number (the 1-based {@code _s##} field) maps directly to the 0-based xlog sfov axis as {@code _s## - 1}. + * Tiles whose SFOV has no (NaN) {@code distance_roi}, or an out-of-range mfov/sfov, are dropped. */ static PointCloud loadPointCloud(final Parameters params, final RenderDataClient renderClient, @@ -457,7 +456,7 @@ static PointCloud loadPointCloud(final Parameters params, params.dataset, offsetX, offsetY); } - // Attach each tile's distance_roi (by mfov + spiral->row-major sfov) and place its center in the voxel frame. + // Attach each tile's distance_roi (by mfov + sfov) and place its center in the voxel frame. final List xs = new ArrayList<>(); final List ys = new ArrayList<>(); final List ds = new ArrayList<>(); @@ -466,14 +465,8 @@ static PointCloud loadPointCloud(final Parameters params, for (final TileBounds tile : tiles) { final String tileId = tile.getTileId(); final int mfov = Integer.parseInt(MultiSemUtilities.getSimpleMfovForTileId(tileId).substring(1)); // m0013 -> 13 - final int spiralSfov = Integer.parseInt(MultiSemUtilities.getSFOVIndexForTileId(tileId)); // _s## spiral, 1-based - // guard the spiral number before the permutation lookup (91-element table) and the mfov before the grid. - if (mfov < 0 || mfov >= nMfov || spiralSfov < 1 || spiralSfov > MultiSemUtilities.NUMBER_OF_TILES_IN_MFOV) { - outOfRange++; - continue; - } - final int sfov = MultiSemUtilities.getRowMajorSFOVIndex(spiralSfov) - 1; // xlog row-major, 0-based - if (sfov < 0 || sfov >= nSfov) { + final int sfov = Integer.parseInt(MultiSemUtilities.getSFOVIndexForTileId(tileId)) - 1; + if (mfov < 0 || mfov >= nMfov || sfov < 0 || sfov >= nSfov) { outOfRange++; continue; } From 9b870e6737e03c04b5999af613d8e1d1d39479b5 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Mon, 13 Jul 2026 13:51:58 -0400 Subject: [PATCH 13/15] Fix mfov numbering --- .../spark/multisem/Wafer6061Inpainter.java | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 10721494a..49f9ea8a6 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -64,9 +64,9 @@ * first layer's tile bounds are fetched from the render web service on the driver, and each tile center (render world * pixels) is mapped to the voxel frame with {@code voxel = center - translate} (the neuroglancer group-level * {@code translate}, i.e. the stack bounding-box min). Each tile carries the {@code distance_roi} of its SFOV, read - * from the xlog for the slab whose {@code id_serial} equals {@code serial}: the tile's mfov (0-based) indexes the xlog - * mfov axis directly, and its sfov number (the 1-based {@code _s##} field) maps directly to the 0-based xlog sfov axis - * as {@code _s## - 1}. + * from the xlog for the slab whose {@code id_serial} equals {@code serial}: the tile's mfov (0-based) maps to the xlog + * mfov axis as {@code render_mfov + 5} (the axis reserves 5 leading always-NaN rows), and its sfov number (the 1-based + * {@code _s##} field) maps to the 0-based xlog sfov axis as {@code _s## - 1}. *

* The client processes the full-resolution ({@code s0}) blocks of the tissue in parallel. For each block it first * applies the cheap ROI-distance filter (interpolating {@code distance_roi} at the block center via inverse-distance @@ -95,6 +95,11 @@ public class Wafer6061Inpainter { // be unique). After slicing the slab axis, distance_roi reduces to a 2-D (sfov, mfov) plane. private static final int DIST_AXIS_SLAB = 2; + // The distance_roi mfov axis reserves 5 leading (always-NaN) positions before the first real mfov: real mfov data + // starts at row 5 for every wafer-60/61 slab (verified constant + contiguous across all 399 finite slabs). Render + // numbers each section's mfovs from 0, so the render mfov maps to the xlog mfov axis position as render_mfov + 5. + private static final int MFOV_ROW_OFFSET = 5; + private static final Logger LOG = LoggerFactory.getLogger(Wafer6061Inpainter.class); public static class Parameters extends CommandLineParameters { @@ -383,8 +388,9 @@ private void runWithSparkContext(final JavaSparkContext sparkContext, * {@code voxel = center - worldToVoxel} ({@code worldToVoxel} is the neuroglancer group-level {@code translate}, * i.e. the stack bounding-box min; if it is null we fall back to the render stack bounds). Each tile carries the * {@code distance_roi} of its SFOV, read from the xlog for the slab whose {@code id_serial} equals - * {@code params.serial}: the tile's mfov (0-based, parsed from the tileId) indexes the xlog mfov axis directly, and - * its sfov number (the 1-based {@code _s##} field) maps directly to the 0-based xlog sfov axis as {@code _s## - 1}. + * {@code params.serial}: the tile's mfov (0-based, parsed from the tileId) maps to the xlog mfov axis as + * {@code render_mfov + 5} (the axis reserves 5 leading always-NaN rows; see {@link #MFOV_ROW_OFFSET}), and its sfov + * number (the 1-based {@code _s##} field) maps to the 0-based xlog sfov axis as {@code _s## - 1}. * Tiles whose SFOV has no (NaN) {@code distance_roi}, or an out-of-range mfov/sfov, are dropped. */ static PointCloud loadPointCloud(final Parameters params, @@ -464,7 +470,8 @@ static PointCloud loadPointCloud(final Parameters params, int outOfRange = 0; for (final TileBounds tile : tiles) { final String tileId = tile.getTileId(); - final int mfov = Integer.parseInt(MultiSemUtilities.getSimpleMfovForTileId(tileId).substring(1)); // m0013 -> 13 + // render mfov is 0-based per section; the xlog mfov axis reserves 5 leading NaN rows (see MFOV_ROW_OFFSET). + final int mfov = Integer.parseInt(MultiSemUtilities.getSimpleMfovForTileId(tileId).substring(1)) + MFOV_ROW_OFFSET; // m0013 -> 13 -> row 18 final int sfov = Integer.parseInt(MultiSemUtilities.getSFOVIndexForTileId(tileId)) - 1; if (mfov < 0 || mfov >= nMfov || sfov < 0 || sfov >= nSfov) { outOfRange++; From e5dfd1f8661de6befc1656aca5e2733c7a3726d0 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Tue, 14 Jul 2026 10:38:16 -0400 Subject: [PATCH 14/15] Inline some variables --- .../spark/multisem/Wafer6061Inpainter.java | 31 +++++++++---------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 49f9ea8a6..238873570 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -403,7 +403,7 @@ static PointCloud loadPointCloud(final Parameters params, final double[][] distBySfovMfov; // [rowMajorSfov0Based][mfov0Based] try (final N5Reader xlog = new N5Factory().openReader(params.xlogPath)) { - final double[] idSerial = read1d(xlog, "id_serial"); + final double[] idSerial = read1d(xlog); final int slabPosition = findSlabPosition(idSerial, params.serial); if (slabPosition < 0) { final long[] available = new long[idSerial.length]; @@ -415,7 +415,7 @@ static PointCloud loadPointCloud(final Parameters params, } final RandomAccessibleInterval distAll = openDoubles(xlog, "distance_roi"); - requireSlabAxis(distAll, DIST_AXIS_SLAB, idSerial.length, "distance_roi"); + requireSlabAxis(distAll, idSerial.length); // distance_roi is scan-independent: slice the slab axis -> 2-D (sfov, mfov). final RandomAccessibleInterval distSlab = Views.hyperSlice(distAll, DIST_AXIS_SLAB, slabPosition); @@ -497,12 +497,12 @@ static PointCloud loadPointCloud(final Parameters params, } /** Fails fast if a hardcoded slab axis does not carry the expected slab count (guards the layout at the top). */ - private static void requireSlabAxis(final RandomAccessibleInterval img, final int axis, final long slabCount, - final String field) { - if (img.numDimensions() <= axis || img.dimension(axis) != slabCount) { - throw new IllegalArgumentException("xlog field '" + field + "' has dimensions " + + private static void requireSlabAxis(final RandomAccessibleInterval img, final long slabCount) { + final int distAxis = Wafer6061Inpainter.DIST_AXIS_SLAB; + if (img.numDimensions() <= distAxis || img.dimension(distAxis) != slabCount) { + throw new IllegalArgumentException("xlog field '" + "distance_roi" + "' has dimensions " + Arrays.toString(img.dimensionsAsLongArray()) + " but the hardcoded layout " + - "expects " + slabCount + " slabs on axis " + axis + + "expects " + slabCount + " slabs on axis " + Wafer6061Inpainter.DIST_AXIS_SLAB + " (see the xlog axis layout at the top of the class)"); } } @@ -517,8 +517,8 @@ private static double[] toArray(final List list) { } /** Reads a 1-D double dataset (blosc-safe: uses the block-not-found overload to avoid the getAttribute NPE). */ - private static double[] read1d(final N5Reader n5, final String dataset) { - final RandomAccessibleInterval img = openDoubles(n5, dataset); + private static double[] read1d(final N5Reader n5) { + final RandomAccessibleInterval img = openDoubles(n5, "id_serial"); final long length = img.dimension(0); final double[] values = new double[(int) length]; final RandomAccess ra = img.randomAccess(); @@ -551,7 +551,7 @@ private static Iterator inpaintPartition(final Iterator bloc LogUtilities.setupExecutorLog4j("inpaint"); final List modified = new ArrayList<>(); - final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(IDW_K); + final KNearestNeighborSearchOnKDTree search = cloud.buildSearch(); try (final N5Reader reader = openN5Reader(params.n5Path); final N5Writer tissueWriter = params.dryRun ? null : openN5Writer(params.n5Path); @@ -574,7 +574,7 @@ private static Iterator inpaintPartition(final Iterator bloc // weighting (the interpolated value is in microns regardless of the voxel-space query units). final double centerX = block.offset[0] + block.dimensions[0] / 2.0; final double centerY = block.offset[1] + block.dimensions[1] / 2.0; - final double roiDistance = cloud.interpolate(search, centerX, centerY, IDW_POWER); + final double roiDistance = cloud.interpolate(search, centerX, centerY); if (! (roiDistance < params.maxRoiDistance)) { logDecision(gridX, gridY, "outside_roi", -1, roiDistance); continue; @@ -975,7 +975,7 @@ int size() { return xs.length; } - KNearestNeighborSearchOnKDTree buildSearch(final int k) { + KNearestNeighborSearchOnKDTree buildSearch() { final List points = new ArrayList<>(xs.length); final List values = new ArrayList<>(xs.length); for (int i = 0; i < xs.length; i++) { @@ -983,14 +983,13 @@ KNearestNeighborSearchOnKDTree buildSearch(final int k) { values.add(new DoubleType(dists[i])); } final KDTree tree = new KDTree<>(values, points); - return new KNearestNeighborSearchOnKDTree<>(tree, Math.min(k, xs.length)); + return new KNearestNeighborSearchOnKDTree<>(tree, Math.min(Wafer6061Inpainter.IDW_K, xs.length)); } /** Inverse-distance-weighted interpolation of the distance value at (x, y). */ double interpolate(final KNearestNeighborSearchOnKDTree search, final double x, - final double y, - final double power) { + final double y) { search.search(new RealPoint(x, y)); final int numNeighbors = search.getK(); double numerator = 0; @@ -1001,7 +1000,7 @@ KNearestNeighborSearchOnKDTree buildSearch(final int k) { if (r == 0.0) { return v; } - final double w = 1.0 / Math.pow(r, power); + final double w = 1.0 / Math.pow(r, Wafer6061Inpainter.IDW_POWER); numerator += w * v; denominator += w; } From 342d6777c2b188acb467a5b425bebf83bbb238c2 Mon Sep 17 00:00:00 2001 From: Michael Innerberger Date: Tue, 14 Jul 2026 10:48:55 -0400 Subject: [PATCH 15/15] Interpolate also partially masked pixels --- .../render/client/spark/multisem/Wafer6061Inpainter.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java index 238873570..938297519 100644 --- a/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java +++ b/render-ws-spark-client/src/main/java/org/janelia/render/client/spark/multisem/Wafer6061Inpainter.java @@ -666,7 +666,7 @@ static InpaintResult inpaintBlock(final RandomAccessibleInterval 0) { + if (maskAccess.setPositionAndGet(local).get() == 255) { value = original; } else { value = zAverage(tissueAccess, local, 0, zMax);