transformation matrix

TransformationMatrix.java

import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.plugin.filter.PlugInFilter;
import ij.process.ImageProcessor;

import java.awt.*;
import java.util.ArrayList;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

public class TransformationMatrix {
ImagePlus im;
	List<int[]> permutations = PermutationHelper.getIndexPermutations(3);
	List<Vertex> points1, points2;
	ImageProcessor cp1, cp2;

	public Matrix correctMatrix = new Matrix(new double[][]{
			{0.013, 1.088, 18.688},
			{-1, -0.05, 127.5}
	});

	@Override
	public int setup(String s, ImagePlus plus)
	{
		this.im = plus;
		return DOES_8G + STACK_REQUIRED;
	}

	@Override
	public void run(ImageProcessor ip)
	{
		int k = im.getStackSize();
		if (k != 2)
		{
			IJ.error("Two images required. Only " + k + " provided.");
			return;
		}
		ImageStack stack = im.getStack();

		ImageProcessor pr1 = stack.getProcessor(1);
		ImageProcessor pr2 = stack.getProcessor(2);

		points1 = IjUtil.collectPoints(pr1);
		points2 = IjUtil.collectPoints(pr2);

		List<Triangle> t1 = getDelaunayTriangulation(points1, pr1.getWidth(), pr1.getHeight());
		List<Triangle> t2 = getDelaunayTriangulation(points2, pr2.getWidth(), pr2.getHeight());

		cp1 = pr1.convertToColorProcessor();
		cp2 = pr2.convertToColorProcessor();

		cp1.setColor(Color.green);
		cp2.setColor(Color.green);

		t1.forEach(t -> IjUtil.draw(t, cp1));
		t2.forEach(t -> IjUtil.draw(t, cp2));

		Matrix m = findTransformationMatrix(t1, t2);
		System.out.println(m);

		//ScoreSet set = getCorrectTransformationMatrix(t1, t2);


		(new ImagePlus("original", cp1)).show();
		(new ImagePlus("transformed", cp2)).show();
	}

	private List<Triangle> getDelaunayTriangulation(List<Vertex> points, int width, int height)
	{
		Triangulation dt = new Triangulation(points, width, height);
		return dt.getDelaunayTriangles();
	}

	private Matrix calculateMatrix(List<Vertex> original, List<Vertex> transformed)
	{
		double x1_1 = original.get(0).coord(0);
		double x1_2 = original.get(1).coord(0);
		double x1_3 = original.get(2).coord(0);

		double y1_1 = original.get(0).coord(1);
		double y1_2 = original.get(1).coord(1);
		double y1_3 = original.get(2).coord(1);

		double x2_1 = transformed.get(0).coord(0);
		double x2_2 = transformed.get(1).coord(0);
		double x2_3 = transformed.get(2).coord(0);

		double y2_1 = transformed.get(0).coord(1);
		double y2_2 = transformed.get(1).coord(1);
		double y2_3 = transformed.get(2).coord(1);

		double d = x1_1 * (y1_3 - y1_2) + x1_2 * (y1_1 - y1_3) + x1_3 * (y1_2 - y1_1);

		double a11 = 1 / d * (y1_1 * (x2_2 - x2_3) + y1_2 * (x2_3 - x2_1) + y1_3 * (x2_1 - x2_2));
		double a12 = 1 / d * (x1_1 * (x2_3 - x2_2) + x1_2 * (x2_1 - x2_3) + x1_3 * (x2_2 - x2_1));
		double a21 = 1 / d * (y1_1 * (y2_2 - y2_3) + y1_2 * (y2_3 - y2_1) + y1_3 * (y2_1 - y2_2));
		double a22 = 1 / d * (x1_1 * (y2_3 - y2_2) + x1_2 * (y2_1 - y2_3) + x1_3 * (y2_2 - y2_1));
		double a13 = 1 / d * (x1_1 * (y1_3 * x2_2 - y1_2 * x2_3)
				+ x1_2 * (y1_1 * x2_3 - y1_3 * x2_1)
				+ x1_3 * (y1_2 * x2_1 - y1_1 * x2_2));
		double a23 = 1 / d * (x1_1 * (y1_3 * y2_2 - y1_2 * y2_3)
				+ x1_2 * (y1_1 * y2_3 - y1_3 * y2_1)
				+ x1_3 * (y1_2 * y2_1 - y1_1 * y2_2));

		return new Matrix(new double[][]{
				{a11, a12, a13},
				{a21, a22, a23}
		});
	}

	private Matrix findTransformationMatrix(List<Triangle> t1, List<Triangle> t2)
	{
		final Matrix[] bestMatrix = {null};
		final Double[] bestDeviation = {null};

		final Triangle[] second = new Triangle[1];
		Triangle original = t1.get(0);
		t2.forEach(transformed -> {

			for (int[] permutation : permutations)
			{
				Matrix trial = calculateMatrix(original.getVertices(), applyPermutation(transformed.getVertices(), permutation));

				final double[] sum = {0};
				points1.stream().map(v -> v.extend(1)).map(trial::multiply).forEach(v -> sum[0] += points2.stream().mapToDouble(v2 -> distanceCheated(v, v2)).min().orElse(0));

				if (bestDeviation[0] == null || sum[0] < bestDeviation[0])
				{
					second[0] = transformed;
					bestDeviation[0] = sum[0];
					bestMatrix[0] = trial;
				}
			}
		});

		// drawing
		cp1.setColor(Color.black);
		cp2.setColor(Color.black);
		IjUtil.draw(original, cp1);
		IjUtil.draw(second[0], cp2);
		cp2.setColor(Color.red);
		Triangle newT = new Triangle(original.getVertices().stream().map(v -> v.extend(1)).map(bestMatrix[0]::multiply).collect(Collectors.toList()));
		System.out.println(newT);
		IjUtil.draw(newT, cp2);
		// draw transformed points
		cp2.setColor(Color.black);
		points1.stream().map(v -> v.extend(1)).map(bestMatrix[0]::multiply).forEach(v -> IjUtil.draw(v, cp2));


		return bestMatrix[0];
	}

	private double distanceCheated(Vertex v1, Vertex v2)
	{
		double a = v2.coord(0) - v1.coord(0);
		double b = v2.coord(1) - v1.coord(1);
		return a * a + b * b;
	}

	private List<Vertex> applyPermutation(ArrayList<Vertex> vertices, int[] perm)
	{
		return IntStream.of(perm).mapToObj(vertices::get).collect(Collectors.toList());
	}
	
}