Flatfield visualization and application

render-app/src/main/java/org/janelia/alignment/destreak/SmoothMaskStreakCorrector.java

@@ -26,7 +26,7 @@
  * The tunable parameters are:
  * <ul>
  * <li>innerCutoff: inner radius of the cutoff profile in px; smaller values may improve de-streaking, values that are too small tinker with the overall intensity of the image</li>
- * <li>bandWidth: width of the gaussian band in px; smaller values may improve de-streaking, values that are too high blur the image in the direction orthogonal to the angle</li>
+ * <li>bandWidth: width of the gaussian band in px; larger values may improve de-streaking, values that are too high blur the image in the direction orthogonal to the angle</li>
  * <li>angle: angle of the band in deg; this should be orthogonal to the streaks (e.g., choose an angle of 0.0 for vertical streaks and 90.0 vor horizontal ones)</li>
  * </ul>
  *

render-ws-java-client/src/main/java/org/janelia/render/client/FlatfieldVisualizationClient.java

@@ -0,0 +1,140 @@
+package org.janelia.render.client;
+
+import org.janelia.alignment.ImageAndMask;
+import org.janelia.alignment.spec.ChannelSpec;
+import org.janelia.alignment.spec.LayoutData;
+import org.janelia.alignment.spec.LeafTransformSpec;
+import org.janelia.alignment.spec.ResolvedTileSpecCollection;
+import org.janelia.alignment.spec.TileSpec;
+import org.janelia.alignment.spec.TransformSpec;
+import org.janelia.alignment.spec.stack.StackMetaData;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import java.awt.Rectangle;
+import java.io.IOException;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.TreeMap;
+
+/**
+ * This is a one-off client for visualizing the estimated flat field of the Hess Lab's wafer_53
+ */
+public class FlatfieldVisualizationClient {
+	// the path to the data (converted to 8-bit PNGs; must be on a shared filesystem accessible to the web server)
+	private static final String DATA_PATH = "/nrs/hess/render/flatfield/flatfields_n2000_8bit";
+	// the format for the files (fill in z and sfov number for the placeholders, in order)
+	private static final String NAME_FORMAT = "flat_field_z%03d_sfov%03d_n2000.png";
+	// the number of z slices
+	private static final int Z_SLICES = 47;
+
+	// stack location
+	private static final String BASE_URL = "http://renderer-dev.int.janelia.org:8080/render-ws/v1";
+	private static final String OWNER = "hess_wafer_53d";
+	private static final String PROJECT = "slab_000_to_009";
+	private static final String STACK = "flatfield_estimate";
+	// reuse metadata from random stack in same project to set up new stack
+	private static final String TEMPLATE_STACK = "s000_m209";
+
+	public static void main(final String[] args) throws IOException {
+		LOG.info("main: setting up target stack {}", STACK);
+		final RenderDataClient dataClient = new RenderDataClient(BASE_URL, OWNER, PROJECT);
+		final StackMetaData projectSpecificMetaData = dataClient.getStackMetaData(TEMPLATE_STACK);
+		dataClient.setupDerivedStack(projectSpecificMetaData, STACK);
+
+		final Map<Integer, TileSpec> tileSpecTemplate = getTemplateLayer(dataClient);
+		for (int z = 1; z <= Z_SLICES; z++) {
+			final ResolvedTileSpecCollection slice = assembleSlice(tileSpecTemplate, z);
+			dataClient.saveResolvedTiles(slice, STACK, (double) z);
+		}
+
+		LOG.info("main: completing target stack {}", STACK);
+		dataClient.setStackState(STACK, StackMetaData.StackState.COMPLETE);
+	}
+
+	private static Map<Integer, TileSpec> getTemplateLayer(final RenderDataClient dataClient) throws IOException {
+		LOG.info("getTemplateLayer: getting template layer from stack {}", TEMPLATE_STACK);
+		final ResolvedTileSpecCollection realTileSpecs = dataClient.getResolvedTiles(TEMPLATE_STACK, 1.0, ".*_000001_.*");
+
+		final Map<Integer, TileSpec> templateTileSpecs = new HashMap<>();
+		final int rowOffset = 10;
+
+		final TileSpec firstTileSpec = realTileSpecs.getTileSpecs().stream().findFirst().orElseThrow();
+		final double height = firstTileSpec.getHeight();
+		final double width = firstTileSpec.getWidth();
+
+		for (final TileSpec tileSpec : realTileSpecs.getTileSpecs()) {
+			final TileSpec templateTileSpec = new TileSpec();
+
+			// set layout
+			final LayoutData originalLayout = tileSpec.getLayout();
+			final LayoutData layoutData = new LayoutData(
+					"1.0",
+					originalLayout.getTemca(),
+					originalLayout.getCamera(),
+					originalLayout.getImageRow() - rowOffset,
+					originalLayout.getImageCol(),
+					originalLayout.getImageCol() * width / 2,
+					(originalLayout.getImageRow() - rowOffset) * height,
+					0.0
+			);
+			templateTileSpec.setLayout(layoutData);
+
+			// set bounds
+			final Rectangle bounds = new Rectangle(templateTileSpec.getLayout().getStageX().intValue(), templateTileSpec.getLayout().getStageY().intValue(), (int) width, (int) height);
+			templateTileSpec.setBoundingBox(bounds, 2000.0);
+			templateTileSpec.setWidth(width);
+			templateTileSpec.setHeight(height);
+
+			// insert at correct position
+			final int sfovNumber = getSfovNumber(tileSpec.getTileId());
+			templateTileSpecs.put(sfovNumber, templateTileSpec);
+		}
+
+		return templateTileSpecs;
+	}
+
+	private static int getSfovNumber(final String tileId) {
+		final String[] parts = tileId.split("_");
+		return Integer.parseInt(parts[2]);
+	}
+
+	private static ResolvedTileSpecCollection assembleSlice(final Map<Integer, TileSpec> templateTileSpecs, final int z) {
+		LOG.info("assembleSlice: assembling slice {}", z);
+		final ResolvedTileSpecCollection slice = new ResolvedTileSpecCollection();
+		for (final Map.Entry<Integer, TileSpec> entry : templateTileSpecs.entrySet()) {
+			final int sfov = entry.getKey();
+			final TileSpec template = entry.getValue();
+
+			// set metadata
+			final TileSpec tileSpec = new TileSpec();
+			tileSpec.setTileId(String.format("z%03d_sfov%03d", z, sfov));
+			tileSpec.setLayout(template.getLayout());
+			tileSpec.setBoundingBox(template.toTileBounds().toRectangle(), 2000.0);
+			tileSpec.setWidth((double) template.getWidth());
+			tileSpec.setHeight((double) template.getHeight());
+			tileSpec.setZ((double) z);
+
+			// transform image
+			final String translationData = String.format("1 0 0 1 %d %d", tileSpec.getLayout().getStageX().intValue(), tileSpec.getLayout().getStageY().intValue());
+			final TransformSpec translation = new LeafTransformSpec("mpicbg.trakem2.transform.AffineModel2D", translationData);
+			tileSpec.addTransformSpecs(List.of(translation));
+
+			// set image
+			final String imagePath = DATA_PATH + "/" + String.format(NAME_FORMAT, z, sfov);
+			final ChannelSpec channelSpec = createChannelSpec(imagePath);
+			tileSpec.addChannel(channelSpec);
+			slice.addTileSpecToCollection(tileSpec);
+		}
+		return slice;
+	}
+
+	private static ChannelSpec createChannelSpec(final String imagePath) {
+		final TreeMap<Integer, ImageAndMask> mipmapLevels = new TreeMap<>();
+		mipmapLevels.put(0, new ImageAndMask(imagePath, null));
+		return new ChannelSpec("default", null, null, mipmapLevels, null, null);
+	}
+
+	private static final Logger LOG = LoggerFactory.getLogger(FlatfieldVisualizationClient.class);
+}

render-ws-java-client/src/main/java/org/janelia/render/client/MultiSemFlatFieldCorrectionClient.java

@@ -0,0 +1,199 @@
+package org.janelia.render.client;
+
+import com.beust.jcommander.Parameter;
+import com.beust.jcommander.ParametersDelegate;
+import ij.IJ;
+import ij.ImagePlus;
+import ij.process.FloatProcessor;
+import ij.process.ImageProcessor;
+import org.janelia.alignment.ImageAndMask;
+import org.janelia.alignment.loader.ImageLoader;
+import org.janelia.alignment.spec.ChannelSpec;
+import org.janelia.alignment.spec.ResolvedTileSpecCollection;
+import org.janelia.alignment.spec.TileSpec;
+import org.janelia.alignment.spec.stack.StackId;
+import org.janelia.alignment.spec.stack.StackMetaData;
+import org.janelia.alignment.util.ImageProcessorCache;
+import org.janelia.render.client.parameter.CommandLineParameters;
+import org.janelia.render.client.parameter.MultiProjectParameters;
+import org.jetbrains.annotations.NotNull;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
+import java.io.IOException;
+import java.nio.file.Path;
+import java.util.List;
+
+/**
+ * Apply a previously computed flat field estimate to correct shading in a Multi-SEM project.
+ * Originally developed for the wafer_53 Multi-SEM dataset.
+ *
+ * @author Michael Innerberger
+ */
+public class MultiSemFlatFieldCorrectionClient {
+
+
+	// make cache large enough to hold all flat field estimates for one layer
+	private final Parameters params;
+	private final ImageProcessorCache cache;
+
+	public static class Parameters extends CommandLineParameters {
+		@ParametersDelegate
+		private final MultiProjectParameters multiProject = new MultiProjectParameters();
+		@Parameter(names = "--flatFieldLocation", description = "Location of the flat field estimates", required = true)
+		private String flatFieldLocation;
+		@Parameter(names = "--outputRoot", description = "Folder to write corrected images to", required = true)
+		private String outputFolder;
+		@Parameter(names = "--flatFieldFormat", description = "Format of the flat field estimates in printf style with placeholders for z-layer and sfov", required = true)
+		private String flatFieldFormat;
+		@Parameter(names = "--targetStackSuffix", description = "Suffix to append to the stack name for the corrected stack", required = true)
+		private String targetStackSuffix = "_corrected";
+		@Parameter(names = "--inputRoot", description = "Root folder for input data; if given, the structure under this root is replicated in the output folder")
+		private String inputRoot = null;
+		@Parameter(names = "--flatFieldConstantFromZ", description = "Maximum z-layer of flat field estimates to consider. All subsequent z-layers get corrected with the maxium, since the estimates can be bad for the last few z-layers. If not given, all z-layers are considered")
+		private Integer flatFieldConstantFromZ = Integer.MAX_VALUE;
+		@Parameter(names = "--cacheSizeGb", description = "Size of the image processor cache in GB (should be enough to hold one layer of flat field estimates)")
+		private double cacheSizeGb = 1.5;
+	}
+
+	public static void main(final String[] args) {
+		final ClientRunner clientRunner = new ClientRunner(args) {
+			@Override
+			public void runClient(final String[] args) throws IOException {
+
+				final Parameters parameters = new Parameters();
+				parameters.parse(args);
+				LOG.info("runClient: entry, parameters={}", parameters);
+
+				final MultiSemFlatFieldCorrectionClient client = new MultiSemFlatFieldCorrectionClient(parameters);
+				client.correctTiles();
+			}
+		};
+		clientRunner.run();
+	}
+
+	public MultiSemFlatFieldCorrectionClient(final Parameters parameters) {
+		this.params = parameters;
+		final double cacheSizeInBytes = 1_000_000_000 * parameters.cacheSizeGb;
+		this.cache = new ImageProcessorCache((long) cacheSizeInBytes, false, false);
+	}
+
+	public void correctTiles() throws IOException {
+		final RenderDataClient renderClient = new RenderDataClient(params.multiProject.baseDataUrl, params.multiProject.owner, params.multiProject.project);
+		final List<StackId> stacks = params.multiProject.stackIdWithZ.getStackIdList(renderClient);
+
+		// TODO: iterate by z-layer first to re-use the flat field estimate for all stacks?
+		for (final StackId stack : stacks) {
+			correctTilesForStack(params.multiProject.baseDataUrl, stack);
+		}
+	}
+
+	public void correctTilesForStack(final String baseDataUrl, final StackId stack) throws IOException {
+			final RenderDataClient stackClient = new RenderDataClient(baseDataUrl, stack.getOwner(), stack.getProject());
+			final List<Double> zValues = stackClient.getStackZValues(stack.getStack());
+			LOG.info("Correcting tiles for {} with {} z-layers", stack, zValues.size());
+
+			final StackMetaData stackMetaData = stackClient.getStackMetaData(stack.getStack());
+			stackClient.setupDerivedStack(stackMetaData, stack.getStack() + params.targetStackSuffix);
+
+			for (final double z : zValues) {
+				final ResolvedTileSpecCollection tileSpecs = stackClient.getResolvedTiles(stack.getStack(), z);
+
+				for (final TileSpec tileSpec : tileSpecs.getTileSpecs()) {
+					final ImageProcessor ip = loadImageTile(tileSpec);
+					final int sfov = extractSfovNumber(tileSpec);
+					final ImageProcessor flatFieldEstimate = loadFlatFieldEstimate(Math.min(z, (double) params.flatFieldConstantFromZ), sfov);
+
+					applyFlatFieldCorrection(ip, flatFieldEstimate);
+
+					final Path newPath = getNewPath(tileSpec);
+					ensureFolderExists(newPath.getParent());
+					patchTileSpec(tileSpec, newPath);
+					saveImage(ip, tileSpec);
+				}
+
+				stackClient.saveResolvedTiles(tileSpecs, stack.getStack() + params.targetStackSuffix, z);
+			}
+			stackClient.setStackState(stack.getStack() + params.targetStackSuffix, StackMetaData.StackState.COMPLETE);
+	}
+
+	private ImageProcessor loadFlatFieldEstimate(final double z, final int sfov) {
+		final Path imagePath = Path.of(params.flatFieldLocation, String.format(params.flatFieldFormat, (int) z, sfov));
+		final String imageUrl = "file:" + imagePath;
+		final ImageLoader.LoaderType loaderType = ImageLoader.LoaderType.IMAGEJ_DEFAULT;
+		return cache.get(imageUrl, 0, false, false, loaderType, null);
+	}
+
+	private int extractSfovNumber(final TileSpec tileSpec) {
+		final String tileId = tileSpec.getTileId();
+		final String[] parts = tileId.split("_");
+		return Integer.parseInt(parts[1].substring(4));
+	}
+
+	/**
+	 * Apply the flat field estimate to the input image.
+	 * @param ip the input image that is altered in place
+	 * @param flatFieldEstimate the flat field estimate to apply
+	 */
+	private void applyFlatFieldCorrection(final ImageProcessor ip, final ImageProcessor flatFieldEstimate) {
+		// convert to 32-bit grayscale (float) for lossless processing
+		final FloatProcessor fp = ip.convertToFloatProcessor();
+
+		for (int i = 0; i < ip.getPixelCount(); i++) {
+			final double a = fp.getf(i);
+			final double b = flatFieldEstimate.getf(i);
+			fp.setf(i, (float) (a / b));
+		}
+
+		// convert back to original bit depth
+		fp.setMinAndMax(0, 255);
+		ip.setPixels(0, fp);
+	}
+
+	private void patchTileSpec(final TileSpec tileSpec, final Path newPath) {
+		final ChannelSpec firstChannel = tileSpec.getAllChannels().get(0);
+		final ImageAndMask originalImage = firstChannel.getFirstMipmapEntry().getValue();
+		final ImageAndMask newImage = originalImage.copyWithImage(newPath.toString(), null, null);
+		firstChannel.putMipmap(0, newImage);
+	}
+
+	private Path getNewPath(final TileSpec tileSpec) {
+		final Path originalPath = Path.of(tileSpec.getTileImageUrl().replaceFirst("file:", ""));
+		if (params.inputRoot != null) {
+			final Path relativePath = Path.of(params.inputRoot).relativize(originalPath);
+			return Path.of(params.outputFolder).resolve(relativePath);
+		} else {
+			return Path.of(params.outputFolder).resolve(originalPath.getFileName());
+		}
+	}
+
+	private void ensureFolderExists(final Path folder) {
+		final boolean folderExists = folder.toFile().exists() || folder.toFile().mkdirs();
+
+		if (!folderExists) {
+			LOG.error("Could not create output folder: {}", folder);
+			System.exit(1);
+		}
+	}
+
+	private ImageProcessor loadImageTile(final TileSpec tileSpec) {
+		final ChannelSpec firstChannelSpec = tileSpec.getAllChannels().get(0);
+		final String tileId = tileSpec.getTileId();
+		final ImageAndMask imageAndMask = firstChannelSpec.getFirstMipmapImageAndMask(tileId);
+
+		return ImageProcessorCache.DISABLED_CACHE.get(imageAndMask.getImageUrl(),
+													  0,
+													  false,
+													  firstChannelSpec.is16Bit(),
+													  imageAndMask.getImageLoaderType(),
+													  imageAndMask.getImageSliceNumber());
+	}
+
+	private void saveImage(final ImageProcessor ip, final TileSpec tileSpec) {
+		final String tileId = tileSpec.getTileId();
+		final ImagePlus imp = new ImagePlus(tileId, ip);
+		IJ.save(imp, tileSpec.getImagePath());
+	}
+
+	private static final Logger LOG = LoggerFactory.getLogger(MultiSemFlatFieldCorrectionClient.class);
+}

render-ws-java-client/src/test/java/org/janelia/render/client/MultiSemFlatFieldCorrectionClientTest.java

@@ -0,0 +1,34 @@
+package org.janelia.render.client;
+
+import org.janelia.render.client.parameter.CommandLineParameters;
+import org.junit.Test;
+
+/**
+ * Tests the {@link MultiSemFlatFieldCorrectionClient} class.
+ *
+ * @author Michael Innerberger
+ */
+public class MultiSemFlatFieldCorrectionClientTest {
+
+    @Test
+    public void testParameterParsing() throws Exception {
+        CommandLineParameters.parseHelp(new MultiSemFlatFieldCorrectionClient.Parameters());
+    }
+
+    public static void main(String[] args) {
+        args = new String[] {
+                "--baseDataUrl", "http://renderer-dev.int.janelia.org:8080/render-ws/v1",
+                "--owner", "hess_wafer_53d",
+                "--project", "slab_070_to_079",
+                "--stackPattern", "s070_m104",
+                "--targetStackSuffix", "_corrected",
+                "--inputRoot", "/nrs/hess/data/hess_wafer_53/raw/imaging/msem",
+                "--outputRoot", "/nrs/hess/render/flatfield/test/corrected",
+                "--flatFieldLocation", "/nrs/hess/render/flatfield/flatfields_n2000",
+                "--flatFieldFormat", "flat_field_z%03d_sfov%03d_n2000.tif",
+        };
+
+        MultiSemFlatFieldCorrectionClient.main(args);
+    }
+
+}