
package net.imagej.ops.fopd.energy.deconvolution;

import net.imagej.ops.OpService;
import net.imagej.ops.fopd.costfunction.CostFunction;
import net.imagej.ops.fopd.operator.FastConvolver;
import net.imagej.ops.fopd.operator.LinearOperator;
import net.imagej.ops.fopd.regularizer.Regularizer;
import net.imagej.ops.fopd.solver.DefaultSolver;
import net.imagej.ops.fopd.solver.Solver;
import net.imagej.ops.fopd.solver.SolverState;
import net.imagej.ops.special.function.AbstractBinaryFunctionOp;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.type.numeric.RealType;
import net.imglib2.view.Views;

import org.scijava.plugin.Parameter;

/**
 * Initialize the {@link LinearOperator} for the deconvolution problem. Actual
 * {@link Solver}s have to implement
 * {@link AbstractDeconvoltuion#getSolverState(RandomAccessibleInterval[])},
 * {@link AbstractDeconvoltuion#getRegularizer(double)} and
 * {@link AbstractDeconvoltuion#getCostFunction(RandomAccessibleInterval[], LinearOperator[], LinearOperator[])}
 * .
 * 
 * @author Tim-Oliver Buchholz, University of Konstanz
 * @param <T>
 */
public abstract class AbstractDeconvoltuion<T extends RealType<T>> extends
	AbstractBinaryFunctionOp<RandomAccessibleInterval<T>[], RandomAccessibleInterval<T>[], RandomAccessibleInterval<T>>
{

	@Parameter
	private int numIt;

	@Parameter
	protected OpService ops;

	public AbstractDeconvoltuion() {
		super();
	}

	@SuppressWarnings({ "unchecked" })
	public RandomAccessibleInterval<T> calculate(
		RandomAccessibleInterval<T>[] input,
		RandomAccessibleInterval<T>[] kernel)
	{

		if (input.length != kernel.length) {
			throw new IllegalArgumentException(
				"Number of input images differs from number of kernels.");
		}

		final Regularizer<T> tgv = getRegularizer(input.length);

		final LinearOperator<T>[] ascentConvolver =
			new LinearOperator[kernel.length];
		final LinearOperator<T>[] descentConvolver =
			new LinearOperator[kernel.length];

		for (int i = 0; i < kernel.length; i++) {
			ascentConvolver[i] = ops.op(FastConvolver.class, input[i],
				kernel[i]);
			descentConvolver[i] = ops.op(FastConvolver.class, input[i], Views
				.invertAxis(Views.invertAxis(ops.copy().rai(kernel[0]), 0), 1));
		}

		final CostFunction<T> cf = getCostFunction(input, ascentConvolver,
			descentConvolver);
		final SolverState<T> state = getSolverState(input);

		final DefaultSolver<T> solver = ops.op(DefaultSolver.class, state, tgv,
			cf, numIt);
		solver.calculate(state);

		return state.getResultImage(0);
	}

	abstract SolverState<T> getSolverState(
		final RandomAccessibleInterval<T>[] input);

	abstract Regularizer<T> getRegularizer(final double numViews);

	abstract CostFunction<T> getCostFunction(
		final RandomAccessibleInterval<T>[] input,
		final LinearOperator<T>[] ascentOperator,
		final LinearOperator<T>[] descentOperator);

}

@@ -2079,16 +2079,58 @@
 el%5Bi%5D);%0A
+%09%09%09if (kernel%5Bi%5D.numDimensions() == 2) %7B%0A%09
 %09%09%09desce
@@ -2194,16 +2194,17 @@
 ews%0A%09%09%09%09
+%09
 .invertA
@@ -2250,21 +2250,235 @@
 nel%5B
-0%5D), 0), 1));
+i%5D), 0), 1));%0A%09%09%09%7D else if (kernel%5Bi%5D.numDimensions() == 3) %7B%0A%09%09%09%09descentConvolver%5Bi%5D = ops.op(FastConvolver.class, input%5Bi%5D,%0A%09%09%09%09%09%09Views.invertAxis(Views.invertAxis(Views.invertAxis(ops.copy().rai(kernel%5Bi%5D), 0), 1), 2));%0A%09%09%09%7D
 %0A%09%09%7D

Fix 3D deconvolution, descent kernel was not correctly flipped.