org.scijava.ops.image.deconvolve.PadAndRichardsonLucy Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of scijava-ops-image Show documentation
Image processing operations for SciJava Ops.
The newest version!
/*
 * #%L
 * Image processing operations for SciJava Ops.
 * %%
 * Copyright (C) 2014 - 2024 SciJava developers.
 * %%
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 * #L%
 */

package org.scijava.ops.image.deconvolve;

import java.util.ArrayList;
import java.util.List;
import java.util.function.BiFunction;

import net.imglib2.Dimensions;
import net.imglib2.FinalDimensions;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.outofbounds.OutOfBoundsConstantValueFactory;
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory;
import net.imglib2.outofbounds.OutOfBoundsFactory;
import net.imglib2.type.NativeType;
import net.imglib2.type.numeric.ComplexType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.util.Util;

import org.scijava.function.Computers;
import org.scijava.function.Functions;
import org.scijava.function.Inplaces;
import org.scijava.ops.spi.Nullable;
import org.scijava.ops.spi.OpDependency;
import org.scijava.progress.Progress;

/**
 * Richardson Lucy function op that operates on (@link RandomAccessibleInterval)
 * (Lucy, L. B. (1974). "An iterative technique for the rectification of
 * observed distributions".)
 *
 * @author Brian Northan
 * @param 
 * @param 
 * @param 
 * @param 
 * @implNote op names='deconvolve.richardsonLucy', priority='100.'
 */
public class PadAndRichardsonLucy & NativeType, O extends RealType & NativeType, K extends RealType & NativeType, C extends ComplexType & NativeType>
	implements
	Functions.Arity9, RandomAccessibleInterval, O, C, Integer, Boolean, Boolean, long[], OutOfBoundsFactory>, RandomAccessibleInterval>
{

	private Computers.Arity1, RandomAccessibleInterval> computeEstimateOp =
		getComputeEstimateOp();

	@OpDependency(name = "deconvolve.normalizationFactor")
	private Inplaces.Arity5_1, Dimensions, Dimensions, RandomAccessibleInterval, RandomAccessibleInterval> normalizer;

	@OpDependency(name = "deconvolve.firstGuess")
	private Functions.Arity3, O, Dimensions, RandomAccessibleInterval> firstGuess;

	@OpDependency(name = "math.multiply")
	private Computers.Arity2, RandomAccessibleInterval, RandomAccessibleInterval> multiplyOp;

	// TODO: can this go in AbstractFFTFilterF?
	@OpDependency(name = "create.img")
	private BiFunction> outputCreator;

	@OpDependency(name = "filter.padInputFFTMethods")
	private Functions.Arity4, Dimensions, Boolean, OutOfBoundsFactory>, RandomAccessibleInterval> padOp;

	@OpDependency(name = "filter.padShiftKernelFFTMethods")
	private BiFunction, Dimensions, RandomAccessibleInterval> padKernelOp;

	@OpDependency(name = "filter.createFFTOutput")
	private Functions.Arity3> createOp;

	@OpDependency(name = "deconvolve.richardsonLucy", hints = {
		"progress.TRACK" })
	private Computers.Arity12, RandomAccessibleInterval, //
			RandomAccessibleInterval, RandomAccessibleInterval, Boolean, //
			Boolean, C, Integer, Boolean, //
			Computers.Arity1, RandomAccessibleInterval>, //
			List>>, //
			RandomAccessibleInterval, RandomAccessibleInterval> richardsonLucyOp;

	private Boolean nonCirculant;

	private Integer maxIterations;

	/**
	 * TODO
	 *
	 * @param input
	 * @param kernel
	 * @param outType
	 * @param complexType
	 * @param maxIterations max number of iterations
	 * @param nonCirculant indicates whether to use non-circulant edge handling
	 * @param accelerate indicates whether or not to use acceleration
	 * @param borderSize
	 * @param obfInput
	 * @return the output
	 */
	@Override
	public RandomAccessibleInterval apply(RandomAccessibleInterval input,
		RandomAccessibleInterval kernel, O outType, C complexType,
		Integer maxIterations, @Nullable Boolean nonCirculant,
		@Nullable Boolean accelerate, @Nullable long[] borderSize,
		@Nullable OutOfBoundsFactory> obfInput)
	{
		// default to circulant
		if (nonCirculant == null) {
			nonCirculant = false;
			this.nonCirculant = nonCirculant;
		}
		else {
			this.nonCirculant = nonCirculant;
		}

		// out of bounds factory will be different depending on if circulant or
		// non-circulant is used
		if (obfInput == null) {
			if (nonCirculant) {
				obfInput = new OutOfBoundsConstantValueFactory<>(Util
					.getTypeFromInterval(input).createVariable());
			}
			else {
				obfInput = new OutOfBoundsMirrorFactory<>(
					OutOfBoundsMirrorFactory.Boundary.SINGLE);
			}
		}

		// default to no acceleration
		if (accelerate == null) accelerate = false;

		this.maxIterations = maxIterations;

		RandomAccessibleInterval output = outputCreator.apply(input, outType);

		final int numDimensions = input.numDimensions();

		// 1. Calculate desired extended size of the image

		final long[] paddedSize = new long[numDimensions];

		if (borderSize == null) {
			// if no border size was passed in, then extend based on kernel size
			for (int d = 0; d < numDimensions; ++d) {
				paddedSize[d] = (int) input.dimension(d) + (int) kernel.dimension(d) -
					1;
			}

		}
		else {
			// if borderSize was passed in
			for (int d = 0; d < numDimensions; ++d) {

				paddedSize[d] = Math.max(kernel.dimension(d) + 2 * borderSize[d], input
					.dimension(d) + 2 * borderSize[d]);
			}
		}

		Progress.defineTotal(0, 1);
		RandomAccessibleInterval paddedInput = padOp.apply(input,
			new FinalDimensions(paddedSize), true, obfInput);

		RandomAccessibleInterval paddedKernel = padKernelOp.apply(kernel,
			new FinalDimensions(paddedSize));

		computeFilter(input, kernel, paddedInput, paddedKernel, output, paddedSize,
			complexType, accelerate);

		return output;
	}

	/**
	 * create a richardson lucy filter
	 */
	private
		Computers.Arity2, RandomAccessibleInterval, RandomAccessibleInterval>
		createFilterComputer(RandomAccessibleInterval unpaddedInput,
			RandomAccessibleInterval unpaddedKernel,
			RandomAccessibleInterval raiExtendedInput,
			RandomAccessibleInterval raiExtendedKernel,
			RandomAccessibleInterval fftImg, RandomAccessibleInterval fftKernel,
			boolean accelerate, RandomAccessibleInterval output)
	{
		C complexType = Util.getTypeFromInterval(fftImg).createVariable();

		// if non-circulant mode, set up the richardson-lucy computer in
		// non-circulant mode and return it
		if (nonCirculant) {
			Inplaces.Arity1> normalizerSimplified = (
				io) -> {
				normalizer.mutate(io, unpaddedInput, unpaddedKernel, fftImg, fftKernel);
			};

			ArrayList>> list =
				new ArrayList<>();

			list.add(normalizerSimplified);

			return (input, kernel, out) -> {
				richardsonLucyOp.compute(input, kernel, fftImg, fftKernel, true, true,
					complexType, maxIterations, accelerate, computeEstimateOp, list,
					firstGuess.apply(raiExtendedInput, Util.getTypeFromInterval(out),
						out), out);
			};
		}

		// return a richardson lucy computer
		return (input, kernel, out) -> {
			richardsonLucyOp.compute(input, kernel, fftImg, fftKernel, true, true,
				complexType, maxIterations, accelerate, computeEstimateOp, null, null,
				out);
		};
	}

	/**
	 * set up and return the compute estimate op. This function can be over-ridden
	 * to implement different types of richardson lucy (like total variation
	 * richardson lucy)
	 *
	 * @return compute estimate op
	 */
	protected
		Computers.Arity1, RandomAccessibleInterval>
		getComputeEstimateOp()
	{
		// TODO: can we now delete RichardsonLucyUpdate?
		return (in, out) -> {
			multiplyOp.compute(out, in, out);
		};
	}

	// TODO: can we move this to AbstractFFTFilterF?
	/**
	 * create FFT memory, create FFT filter and run it
	 */
	private void computeFilter(final RandomAccessibleInterval input,
		final RandomAccessibleInterval kernel,
		final RandomAccessibleInterval paddedInput,
		final RandomAccessibleInterval paddedKernel,
		RandomAccessibleInterval output, long[] paddedSize, C complexType,
		boolean accelerate)
	{

		RandomAccessibleInterval fftInput = createOp.apply(new FinalDimensions(
			paddedSize), complexType, true);

		RandomAccessibleInterval fftKernel = createOp.apply(new FinalDimensions(
			paddedSize), complexType, true);

		// TODO: in this case it is difficult to match the filter op in the
		// 'initialize' as we don't know the size yet, thus we can't create
		// memory
		// for the FFTs
		Computers.Arity2, RandomAccessibleInterval, RandomAccessibleInterval> filter =
			createFilterComputer(input, kernel, paddedInput, paddedKernel, fftInput,
				fftKernel, accelerate, output);

		filter.compute(paddedInput, paddedKernel, output);
	}

}