/*====================================================================	
| Version: July 29, 2015
\===================================================================*/

/*====================================================================
Minimum spanning tree based alignment of 3D image sequence

Nilanjan Ray
Computing Science
University of Alberta
Canada
\===================================================================*/

/*====================================================================

Description:

This program aligns a sequences of 3D images by minimum spanning tree-based image alignment algorithm MISTICA; If you use our program in your work, 
please cite our relevant papers. To use this program you need to have a sequence of images open as a hyperstack in ImageJ. The program will modify 
the images in the hyperstack.

User input parameters:

Transformation: Choose "Translation" or "Rigid Body" to get best coarse alignment results.
Graph width: If "Exact" option is used for the Method, then use a smaller value here (say, between 1 to 5)
Anchor selection: Default "Automatic" option obtains the best result for most cases; For the "First Z-stack" option, the 
					current slice in the stack will be treated as the anchor image to which other images are registered.
Method: "Approximate" option is fast. "Exact" option is much slower. Exact method might lead to more accurate results.
Add borders: Whether the user wants to increase the image size after transformations. 

Dependencies:

This plugin requires these two packages:
(1) turboreg (http://bigwww.epfl.ch/thevenaz/turboreg/) and 
(2) jgrapht (http://jgrapht.org/)

Acknowledgements: 

The plugin uses some methods taken from the StackReg plugin (http://bigwww.epfl.ch/thevenaz/stackreg/)


\===================================================================*/

// ImageJ
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.WindowManager;
import ij.gui.GUI;
import ij.gui.ImageWindow;
import ij.gui.GenericDialog;
import ij.io.FileSaver;
import ij.plugin.PlugIn;
import ij.process.*;

// Java 1.1
import java.awt.BorderLayout;
import java.awt.Button;
import java.awt.Dialog;
import java.awt.FlowLayout;
import java.awt.Frame;
import java.awt.Insets;
import java.awt.Label;
import java.awt.Panel;
import java.awt.TextArea;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.awt.image.IndexColorModel;
import java.lang.reflect.InvocationTargetException;
import java.lang.reflect.Method;


// jGraphT
import java.net.*;
import java.util.*;

import org.jgrapht.*;
import org.jgrapht.graph.*;
import org.jgrapht.alg.*;
import org.jgrapht.traverse.*;

// For visualization
/*
import javax.swing.JFrame;
import org.jgraph.JGraph;
import org.jgraph.graph.DefaultGraphCell;
import org.jgraph.graph.GraphConstants;
import org.jgrapht.ext.JGraphModelAdapter;
*/
/*====================================================================
|	MSTHyperStackReg3_
\===================================================================*/
public class MSTHyperStackReg3_
	implements
		PlugIn

{ /* class MSTHyperStackReg3_ */

/*....................................................................
	private variables
....................................................................*/
private static final double TINY =
	(double)Float.intBitsToFloat((int)0x33FFFFFF);

private int graphWidth = 5; // width of graph
private double[][][][] transformationMatrix = null; 
SimpleWeightedGraph<Integer, DefaultWeightedEdge>  graph = null;
private int width=0; // image width
private int height=0; // image height
private int nSlice=0; // number of slices in each time point
private int nFrame=0; // total number of time points

/*....................................................................
	PlugIn methods
....................................................................*/
/*------------------------------------------------------------------*/
public void run (
	final String arg
) {
	// prelims and dialog box
	Runtime.getRuntime().gc();
	final ImagePlus imp = WindowManager.getCurrentImage();
	if (imp == null) {
		IJ.error("No image available");
		return;
	}
	if (imp.getStack().isRGB() || imp.getStack().isHSB()) {
		IJ.error("Unable to process either RGB or HSB stacks");
		return;
	}
	GenericDialog gd = new GenericDialog("MSTHyperStackReg3_");
	final String[] transformationItem = {
		"Translation",
		"Rigid Body",
		"Scaled Rotation",
		"Affine"
	};
	final String[] anchor = {
			"Automatic",
			"First Z-stack"
	};
	final String[] method = {
			"Approximate",
			"Exact"
	};
	final String[] border = {
			"No",
			"Yes"
	};
	gd.addChoice("Transformation:", transformationItem, "Rigid Body");
	gd.addNumericField("Intra-Stack graph width:", 0, 0); // default value is 0, meaning that intra-stack registration can be skipped.
	gd.addNumericField("Inter stack graph width:", graphWidth, 0);
	gd.addChoice("Anchor selection", anchor, "Automatic");
	gd.addChoice("Method", method, "Approximate");
	gd.addChoice("Add border?", border, "No");

	gd.showDialog();
	if (gd.wasCanceled()) {
		return;
	}
	final int transformation = gd.getNextChoiceIndex();
	final int intraStackGraphWidth = (int)gd.getNextNumber();
	final int interStackGraphWidth = (int)gd.getNextNumber();
	final int anchorChoice = gd.getNextChoiceIndex();
	final int methodChoice = gd.getNextChoiceIndex();
	final int borderChoice = gd.getNextChoiceIndex();
	
	// get dimensions of the stack / hyperstack and set variables
	int[] dim = imp.getDimensions();		
	if(dim[2] != 1){
		IJ.showMessage("MSTHyperStackReg3_","Requires # channels to be 1");
		return;
	}
	width=dim[0];
	height=dim[1];
	nSlice=dim[3];
	nFrame=dim[4];

	int targetSlice = imp.getCurrentSlice();
	double[][] anchorPoints = computeAnchorPoints(transformation,width,height);
	// color weights
	ImageStack stack = imp.getStack();
	double[] colorWeights = null;
	switch (imp.getType()) {
		case ImagePlus.COLOR_256:
		case ImagePlus.COLOR_RGB: {
			colorWeights = getColorWeightsFromPrincipalComponents(stack,targetSlice);
		}
	}

	if (intraStackGraphWidth>=1){ // skip intra-stack registration if graph width is zero or less
		// STEP 1: Intra-stack registration: stacks are independently registered
		// graphWidth =1 for intra stack registration
		graphWidth = intraStackGraphWidth;
		// create transformation matrix to hold registration information
		transformationMatrix = new double[nSlice*nFrame][2*graphWidth][3][3];
		// bounding box
		int maxt=0; int maxb=height; int maxl=0; int maxr=width;
		// Lists for each stack
		List visitedStack = new ArrayList<List>();
		List parentsStack = new ArrayList<List>();
		List transformMatStack = new ArrayList<List>();

		for(int nf=1; nf<=nFrame; nf++){ // for each stack
	
			// create a graph with vertices as the images in the stack
			graph = new SimpleWeightedGraph<Integer, DefaultWeightedEdge>(DefaultWeightedEdge.class); 
			for(int i=1; i<=nSlice;i++){
				graph.addVertex((nf-1)*nSlice+i);
			}
			// compute transformation, registration error between image pairs and set graph edge weights
			for (int i=1; i<=nSlice;i++){
				ImagePlus target = computeTarget(imp, stack, (nf-1)*nSlice+i, colorWeights, width, height);
				if(i+1<=nSlice){
					if(methodChoice==0){
						// approximate graph creation; fast and may be inaccurate
						appoximateGraphCreate(imp, stack, (nf-1)*nSlice+i, (nf-1)*nSlice+i+1, colorWeights, width, height);
					}else{
					// exact graph creation; slow, but accurate
						registerSlice(target, imp, stack, width, height,transformation, anchorPoints, colorWeights, (nf-1)*nSlice+i, (nf-1)*nSlice+i+1, true);
					}
				}
				IJ.showStatus("Intra-stack registration. Building graph for z stack "+ IJ.d2s(nf,0) +" : " + IJ.d2s(i,0) + " out of " + IJ.d2s(nSlice,0) + " images processed.");
			}
	
			IJ.showStatus("Computing minimum spanning tree...");
			// Compute minimum spanning tree
			KruskalMinimumSpanningTree<String, DefaultWeightedEdge>  kmsTree = new KruskalMinimumSpanningTree(graph);
			// construct the minimum weighted spanning tree returned by KruskalMinimumSpanningTree
			UndirectedSubgraph<Integer, DefaultWeightedEdge>  tree = new UndirectedSubgraph(graph, graph.vertexSet(), kmsTree.getMinimumSpanningTreeEdgeSet());
			IJ.showStatus("Computing minimum spanning tree...done.");

	// for visualization of the tree
    //	JFrame frame = new JFrame();
    //	frame.setSize(400, 400);
    //	JGraph jgraph = new JGraph(new JGraphModelAdapter(tree));
    //	frame.getContentPane().add(jgraph);
    //	frame.setVisible(true);

	
			// Always automatic anchor image computation for intra-stack registration - it should be half way on the longest path
			// approximately, the longest path is between the first and the last node
			IJ.showStatus("Finding out anchor image...");
			DijkstraShortestPath<Integer, DefaultWeightedEdge> spath = new DijkstraShortestPath<Integer, DefaultWeightedEdge>(tree,(nf-1)*nSlice+1,nf*nSlice);
			List list = Graphs.getPathVertexList(spath.getPath());
			int mid = list.size();
			mid = mid/2;
			mid = (Integer)list.get(mid);
			IJ.showStatus("Finding out anchor image...done.");
			targetSlice = mid;
			imp.setSlice(targetSlice);
	
			// breadth first tree traversal and some book keeping
			BreadthFirstIterator<Integer, DefaultWeightedEdge> bfs = new BreadthFirstIterator<Integer, DefaultWeightedEdge>(tree, targetSlice);
			List visited = new ArrayList<Integer>();
			List parents = new ArrayList<Integer>();
			List transformMat = new ArrayList<Double[][]>();
			int next = bfs.next();
			visited.add(next);
			parents.add(0);
			double[][] tmtx = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{0.0, 0.0, 1.0}};
			transformMat.add(tmtx);
			while(bfs.hasNext()){
				next = bfs.next();
				visited.add(next);
				List neighbors = Graphs.neighborListOf(tree, next);
				Object[] neighborsObj = neighbors.toArray();
				int parent=0;
				for(int i=0; i<neighbors.size(); i++){
					int j = (Integer)neighborsObj[i];
					if(visited.contains(j)){
						parent=j;
						break;
					}
				}
				parents.add(parent);
				int k = parent;
				int l = next;
				int m=-1;
				if(k<l)
					m = l-k-1;
				else
					m = graphWidth + k-l-1;
					double[][] tmpmatrix = {
						{transformationMatrix[k-1][m][0][0], transformationMatrix[k-1][m][0][1], transformationMatrix[k-1][m][0][2]},
						{transformationMatrix[k-1][m][1][0], transformationMatrix[k-1][m][1][1], transformationMatrix[k-1][m][1][2]},
						{transformationMatrix[k-1][m][2][0], transformationMatrix[k-1][m][2][1], transformationMatrix[k-1][m][2][2]}
				};
				double[][] tmatrix = (double[][])transformMat.get(visited.indexOf(parent));
				tmatrix = matrixMultiply(tmpmatrix,tmatrix);
				transformMat.add(tmatrix);
			}
	
			// compute bounding box for the image stack
			for(int i=1; i<visited.size(); i++){
				double[][] tmatrix = (double[][])transformMat.get(i);
				int[] bbox = computeBoundingBox(tmatrix, width, height, maxt, maxb, maxr, maxl);
				maxt = bbox[0]; maxb = bbox[1]; maxr = bbox[2]; maxl = bbox[3];
			}
		
			// keep the lists 
			visitedStack.add(visited);
			parentsStack.add(parents);
			transformMatStack.add(transformMat);

		}
	
		if (borderChoice==1){
			// resizing canvas
			maxt +=10; maxb +=10; maxr +=10; maxl +=10; // add some extra pixels
			// zeropad images in the stack
			if(maxr>width)
				width = maxr;
			if(maxb>height)
				height = maxb;
			IJ.run("Canvas Size...", "width="+width+" height="+height+" position=Top-Left zero");
			width += maxl;
			height += maxt;
			IJ.run("Canvas Size...", "width="+width+" height="+height+" position=Bottom-Right zero");
	
			// recompute anchor points
			anchorPoints = computeAnchorPoints(transformation,width,height);
		}
	
		for(int nf=1; nf<=nFrame; nf++){
			List visited = (List)visitedStack.get(nf-1);
			List parents = (List)parentsStack.get(nf-1);
			List transformMat = (List)transformMatStack.get(nf-1);
			// transform the images
			for(int i=1; i<visited.size(); i++){
				// transform the slice first according to already computed transformation
				int v = (Integer)visited.get(i);
				double[][] tmatrix = (double[][])transformMat.get(i);
				ImagePlus source = transformSlice(imp, width, height,transformation, tmatrix, anchorPoints, v);		
				IJ.showStatus("Intra-stack registration. Registering and updating z stack "+ IJ.d2s(nf,0) +" : " + IJ.d2s(i,0) + " out of " + IJ.d2s(nSlice,0) + " images processed.");
			}
			imp.setSlice(targetSlice);
	
			// finetune registration
			stack = imp.getStack();
			for(int i=1; i<visited.size(); i++){
				int t = (Integer)parents.get(i);
				int s = (Integer)visited.get(i);
				ImagePlus target = computeTarget(imp, stack, t, colorWeights, width, height);
				registerSlice(target, imp, stack, width, height,transformation, anchorPoints, colorWeights, t, s, false);
				IJ.showStatus("Intra-stack registration. Finetuning registration for z stack "+ IJ.d2s(nf,0) +" : " + IJ.d2s(i+1,0) + " out of " + IJ.d2s(nSlice,0) + " images processed.");
			}	
			for(int i=1; i<visited.size(); i++){
				int t = (Integer)parents.get(i);
				int s = (Integer)visited.get(i);
				int m=-1;
				if(t<s)
					m = s-t-1;
				else
					m = graphWidth + t-s-1;
				double[][] tmpmatrix = {
					{transformationMatrix[t-1][m][0][0], transformationMatrix[t-1][m][0][1], transformationMatrix[t-1][m][0][2]},
					{transformationMatrix[t-1][m][1][0], transformationMatrix[t-1][m][1][1], transformationMatrix[t-1][m][1][2]},
					{transformationMatrix[t-1][m][2][0], transformationMatrix[t-1][m][2][1], transformationMatrix[t-1][m][2][2]}
				};
				double[][] tmatrix = (double[][])transformMat.get(visited.indexOf(t));
				tmatrix = matrixMultiply(tmpmatrix,tmatrix);
				ImagePlus source = transformSlice(imp, width, height,transformation, tmatrix, anchorPoints, s);
				tmatrix = (double[][])transformMat.set(i,tmatrix);
				IJ.showStatus("Intra-stack registration. Updating z stack "+ IJ.d2s(nf,0) +" : " + IJ.d2s(i+1,0) + " out of " + IJ.d2s(nSlice,0) + " images processed.");
			}
			imp.setSlice(targetSlice);
		}
	} 

// STEP 2: Course Inter stack registration: only a global transformation is computed for an entire stack
	graphWidth = interStackGraphWidth;
	IJ.run("Z Project...", "start=1 stop="+ nSlice +" projection=[Max Intensity] all");
	final ImagePlus imp2 = WindowManager.getCurrentImage();
	
	int[] dim2 = imp2.getDimensions();		
	targetSlice = imp2.getCurrentSlice();
	ImageStack stack2 = imp2.getStack();
	// create transformation matrix to hold registration information
	transformationMatrix = new double[imp2.getStackSize()][2*graphWidth][3][3];
	// create a graph with vertices as the images in the stack
	graph = new SimpleWeightedGraph<Integer, DefaultWeightedEdge>(DefaultWeightedEdge.class); 
	for(int i=1; i<=imp2.getStackSize();i++){
		graph.addVertex(i);
	}
    // colorweights
	double[] colorWeights2 = null;
	switch (imp2.getType()) {
		case ImagePlus.COLOR_256:
		case ImagePlus.COLOR_RGB: {
			colorWeights2 = getColorWeightsFromPrincipalComponents(stack,targetSlice);
		}
	}
	ImageStack newStack = null;
	if(methodChoice==0){
		// preprocess stack: most expensive component when blurring is applied
		newStack = preprocess(imp2,stack2);
	}
	for (int i=1; i<=imp2.getStackSize();i++){
		ImagePlus target = computeTarget(imp2, stack2, i, colorWeights2, width, height);
			//imp.setSliceWithoutUpdate(i);
		for (int j=1; j<=graphWidth; j++){
			if((i+j)<=imp2.getStackSize()){
				if(methodChoice==0){
					// approximate graph creation; fast and may be inaccurate
					appoximateGraphCreate(imp2, newStack, i, i+j, colorWeights2, newStack.getWidth(), newStack.getHeight());
				}else{
					// exact graph creation; slow, but accurate
					registerSlice(target, imp2, stack2, width, height,transformation, anchorPoints, colorWeights2, i, i+j, true);
				}
			}
		}
		IJ.showStatus("Inter stack registration. Building graph: " + IJ.d2s(i,0) + " out of " + IJ.d2s(imp2.getStackSize(),0) + " images processed.");
	}
	IJ.showStatus("Computing minimum spanning tree...");
	// Compute minimum spanning tree
	KruskalMinimumSpanningTree<String, DefaultWeightedEdge>  kmsTree = new KruskalMinimumSpanningTree(graph);
	// construct the minimum weighted spanning tree returned by KruskalMinimumSpanningTree
	UndirectedSubgraph<Integer, DefaultWeightedEdge>  tree = new UndirectedSubgraph(graph, graph.vertexSet(), kmsTree.getMinimumSpanningTreeEdgeSet());
	IJ.showStatus("Computing minimum spanning tree...done.");
	
	targetSlice = 1;
	if(anchorChoice==0){// if anchor choice is automatic
		IJ.showStatus("Finding out anchor image...");
		DijkstraShortestPath<Integer, DefaultWeightedEdge> spath = new DijkstraShortestPath<Integer, DefaultWeightedEdge>(tree,1,imp2.getStackSize());
		List list = Graphs.getPathVertexList(spath.getPath());
		int mid = list.size();
		mid = mid/2;
		mid = (Integer)list.get(mid);
		targetSlice = mid;
		IJ.showStatus("Finding out anchor image...done.");
	}
	imp2.setSlice(targetSlice);

	// breadth first tree traversal and some book keeping
	BreadthFirstIterator<Integer, DefaultWeightedEdge> bfs = new BreadthFirstIterator<Integer, DefaultWeightedEdge>(tree, targetSlice);
	List visited = new ArrayList<Integer>();
	List parents = new ArrayList<Integer>();
	List transformMat = new ArrayList<Double[][]>();
	int next = bfs.next();
	visited.add(next);
	parents.add(0);
	double[][] tmtx = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{0.0, 0.0, 1.0}};
	transformMat.add(tmtx);
	while(bfs.hasNext()){
		next = bfs.next();
		visited.add(next);
		List neighbors = Graphs.neighborListOf(tree, next);
		Object[] neighborsObj = neighbors.toArray();
		int parent=0;
		for(int i=0; i<neighbors.size(); i++){
			int j = (Integer)neighborsObj[i];
			if(visited.contains(j)){
				parent=j;
				break;
			}
		}
		parents.add(parent);
		int k = parent;
		int l = next;
		if(methodChoice==0){ // approximate graph
			ImagePlus target = computeTarget(imp2, stack2, k, colorWeights2, width, height);
			registerSlice(target, imp2, stack2, width, height,transformation, anchorPoints, colorWeights2, k, l, false);
		}
		int m=-1;
		if(k<l)
			m = l-k-1;
		else
			m = graphWidth + k-l-1;
		double[][] tmpmatrix = {
				{transformationMatrix[k-1][m][0][0], transformationMatrix[k-1][m][0][1], transformationMatrix[k-1][m][0][2]},
				{transformationMatrix[k-1][m][1][0], transformationMatrix[k-1][m][1][1], transformationMatrix[k-1][m][1][2]},
				{transformationMatrix[k-1][m][2][0], transformationMatrix[k-1][m][2][1], transformationMatrix[k-1][m][2][2]}
			};
		double[][] tmatrix = (double[][])transformMat.get(visited.indexOf(parent));
		tmatrix = matrixMultiply(tmpmatrix,tmatrix);
		transformMat.add(tmatrix);
	}
	// transform the images
	for(int i=1; i<visited.size(); i++){
		// transform the slice first according to already computed transformation
		int v = (Integer)visited.get(i);
		double[][] tmatrix = (double[][])transformMat.get(i);
		ImagePlus source = transformSlice(imp2, width, height,transformation, tmatrix, anchorPoints, v);
		// transform images in the original stack
		for(int ns=1; ns<=nSlice; ns++){
			source = transformSlice(imp, width, height,transformation, tmatrix, anchorPoints, (v-1)*nSlice+ns);			
		}
		IJ.showStatus("Inter stack registration. Registering and updating stack: " + IJ.d2s(i,0) + " out of " + IJ.d2s(imp2.getStackSize(),0) + " stacks processed.");
	}
	imp2.setSlice(targetSlice);
	// finetune registration
	stack2 = imp2.getStack();
	for(int i=1; i<visited.size(); i++){
		int t = (Integer)parents.get(i);
		int s = (Integer)visited.get(i);
		ImagePlus target = computeTarget(imp2, stack2, t, colorWeights2, width, height);
		registerSlice(target, imp2, stack2, width, height,transformation, anchorPoints, colorWeights2, t, s, false);
		IJ.showStatus("Inter stack registration. Finetuning registration: " + IJ.d2s(i+1,0) + " out of " + IJ.d2s(imp2.getStackSize(),0) + " images processed.");
	}
	for(int i=1; i<visited.size(); i++){
		int t = (Integer)parents.get(i);
		int s = (Integer)visited.get(i);
		int m=-1;
		if(t<s)
			m = s-t-1;
		else
			m = graphWidth + t-s-1;
		double[][] tmpmatrix = {
				{transformationMatrix[t-1][m][0][0], transformationMatrix[t-1][m][0][1], transformationMatrix[t-1][m][0][2]},
				{transformationMatrix[t-1][m][1][0], transformationMatrix[t-1][m][1][1], transformationMatrix[t-1][m][1][2]},
				{transformationMatrix[t-1][m][2][0], transformationMatrix[t-1][m][2][1], transformationMatrix[t-1][m][2][2]}
			};
		double[][] tmatrix = (double[][])transformMat.get(visited.indexOf(t));
		tmatrix = matrixMultiply(tmpmatrix,tmatrix);
		ImagePlus source = transformSlice(imp2, width, height,transformation, tmatrix, anchorPoints, s);
		// transform images in the original stack
		for(int ns=1; ns<=nSlice; ns++){
			source = transformSlice(imp, width, height,transformation, tmatrix, anchorPoints, (s-1)*nSlice+ns);			
		}
		tmatrix = (double[][])transformMat.set(i,tmatrix);
		IJ.showStatus("Inter stack registration. Updating stack: " + IJ.d2s(i+1,0) + " out of " + IJ.d2s(imp2.getStackSize(),0) + " stacks processed.");
	}
	imp2.setSlice(targetSlice);
		
	
	IJ.showStatus("MSTHyperStackReg is done.");
} /* end run */

/*....................................................................
	private methods
....................................................................*/

/*------------------------------------------------------------------*/

private ImageStack preprocess(final ImagePlus imp, ImageStack stack){
	// preprocess stack: create a new stack of 1/resizeFactor the width and the height of the input stack
	final int resizeFactor = 4;
	final double sigmaGauss = 2; 
	ImageStack stackOut = new ImageStack((int)(stack.getWidth()/resizeFactor),(int)(stack.getHeight()/resizeFactor));
	for(int n=1; n<=imp.getStackSize(); n++){
		ImageProcessor im = stack.getProcessor(n);
		ImageProcessor imDuplicate = im.duplicate();
		if(sigmaGauss>=0.5)
			imDuplicate.blurGaussian(sigmaGauss); // Gaussian smoothing filter to suppress noise
		imDuplicate.setInterpolationMethod(ImageProcessor.BILINEAR);
		ImageProcessor im1 = imDuplicate.resize((int)(stack.getWidth()/resizeFactor),(int)(stack.getHeight()/resizeFactor));
		stackOut.addSlice(im1);
		IJ.showProgress(n, imp.getStackSize());
		IJ.showStatus("Preprocessing: " + IJ.d2s(n,0) + " out of " + IJ.d2s(imp.getStackSize(),0) + " images processed.");
	}
	return stackOut;
}

/*------------------------------------------------------------------*/
private int[] computeBoundingBox(final double[][] tmatrix, final int width, final int height, int maxt, int maxb, int maxr, int maxl){
	int[] bbox = new int[4];
	// corner 1
	double x = 0*tmatrix[0][0] - 0*tmatrix[0][1] + tmatrix[0][2];
	double y = 0*tmatrix[1][0] - 0*tmatrix[1][1] + tmatrix[1][2];
	y = -y;
	if(x>maxr){
		maxr = (int)x;
	}else{
		if(-x>maxl)
			maxl = -(int)x;
	}
	if(y>maxb){
		maxb = (int)y;
	}else{
		if(-y>maxt)
			maxt = -(int)y;
	}
	// corner 2
	x = (width-1)*tmatrix[0][0] - 0*tmatrix[0][1] + tmatrix[0][2];
	y = (width-1)*tmatrix[1][0] - 0*tmatrix[1][1] + tmatrix[1][2];
	y = -y;
	if(x>maxr){
		maxr = (int)x;
	}else{
		if(-x>maxl)
			maxl = -(int)x;
	}
	if(y>maxb){
		maxb = (int)y;
	}else{
		if(-y>maxt)
			maxt = -(int)y;
	}
	// corner 3
	x = (width-1)*tmatrix[0][0] - (height-1)*tmatrix[0][1] + tmatrix[0][2];
	y = (width-1)*tmatrix[1][0] - (height-1)*tmatrix[1][1] + tmatrix[1][2];
	y = -y;
	if(x>maxr){
		maxr = (int)x;
	}else{
		if(-x>maxl)
			maxl = -(int)x;
	}
	if(y>maxb){
		maxb = (int)y;
	}else{
		if(-y>maxt)
			maxt = -(int)y;
	}
	// corner 4
	x = (0)*tmatrix[0][0] - (height-1)*tmatrix[0][1] + tmatrix[0][2];
	y = (0)*tmatrix[1][0] - (height-1)*tmatrix[1][1] + tmatrix[1][2];
	y = -y;
	if(x>maxr){
		maxr = (int)x;
	}else{
		if(-x>maxl)
			maxl = -(int)x;
	}
	if(y>maxb){
		maxb = (int)y;
	}else{
		if(-y>maxt)
			maxt = -(int)y;
	}
	bbox[0] = maxt; bbox[1] = maxb; bbox[2] = maxr; bbox[3] = maxl;
	return bbox;
}


/*------------------------------------------------------------------*/
private double mad(ImagePlus im1, ImagePlus im2){
	
	double pixelMad=0.0;
	switch (im1.getType()) {
		case ImagePlus.GRAY8: {
			int len = im1.getWidth()*im1.getHeight();
			byte[] p1 = (byte[])im1.getProcessor().getPixels();
			float[] p2 = (float[])im2.getProcessor().getPixels();
			for(int i=0; i<len; i++){
				float diff = (float)p1[i]-p2[i];
				if(diff<0){ diff = -diff;}
				pixelMad += (double)diff/(double)len;
			}
			return pixelMad;
		}
		case ImagePlus.GRAY16: {
			int len = im1.getWidth()*im1.getHeight();
			short[] p1 = (short[])im1.getProcessor().getPixels();
			float[] p2 = (float[])im2.getProcessor().getPixels();
			for(int i=0; i<len; i++){
				float diff = (float)p1[i]-p2[i];
				if(diff<0){ diff = -diff;}
				pixelMad += (double)diff/(double)len;
			}
			return pixelMad;
		}
		case ImagePlus.GRAY32: {
			int len = im1.getWidth()*im1.getHeight();
			float[] p1 = (float[])im1.getProcessor().getPixels();
			float[] p2 = (float[])im2.getProcessor().getPixels();
			for(int i=0; i<len; i++){
				float diff = p1[i]-p2[i];
				if(diff<0){ diff = -diff;}
				pixelMad += (double)diff/(double)len;
			}
			return pixelMad;
		}
		default: {
			IJ.error("Computing cost: Unexpected image type");
			return -1;
		}
	} // end switch

}


/*------------------------------------------------------------------*/
//Approximating graph creation by resized image difference

private void appoximateGraphCreate(ImagePlus imp, ImageStack stack, int t, int s, double[] colorWeights, int width, int height){
	// create edge and add cost to the graph
	
	ImagePlus im1 = computeTarget(imp, stack, t, colorWeights, width, height);
	ImagePlus im2 = computeTarget(imp, stack, s, colorWeights, width, height);
	double cost = mad(im1, im2);
	if(t<s){
		DefaultWeightedEdge e1 = graph.addEdge(t, s);
		graph.setEdgeWeight(e1, cost); 
	}
	else{
		DefaultWeightedEdge e1 = graph.addEdge(s, t);
		graph.setEdgeWeight(e1, cost); 
	}
}

/*-------------------------------------------------------------------*/
private ImagePlus computeTarget(ImagePlus imp, ImageStack stack, int i, double[] colorWeights, int width, int height){
	ImagePlus target = null;
	//imp.setSliceWithoutUpdate(i);
	switch (imp.getType()) {
		case ImagePlus.COLOR_256:
		case ImagePlus.COLOR_RGB: {
			target = getGray32("StackRegTarget", stack, i, colorWeights);
			break;
		}
		case ImagePlus.GRAY8: {
			target = new ImagePlus("StackRegTarget",
				new ByteProcessor(width, height, new byte[width * height],
				stack.getProcessor(i).getColorModel()));
			target.getProcessor().copyBits(
				stack.getProcessor(i), 0, 0, Blitter.COPY);
			break;
		}
		case ImagePlus.GRAY16: {
			target = new ImagePlus("StackRegTarget",
				new ShortProcessor(width, height, new short[width * height],
				stack.getProcessor(i).getColorModel()));
			target.getProcessor().copyBits(
				stack.getProcessor(i), 0, 0, Blitter.COPY);
			break;
		}
		case ImagePlus.GRAY32: {
			target = new ImagePlus("StackRegTarget",
				new FloatProcessor(width, height, new float[width * height],
				stack.getProcessor(i).getColorModel()));
			target.getProcessor().copyBits(
				stack.getProcessor(i), 0, 0, Blitter.COPY);
			break;
		}
		default: {
			IJ.error("Unexpected image type");
			return null;
		}
	}
	return target;
}
/*-------------------------------------------------------------------*/
private double[][] computeAnchorPoints(final int transformation, final int width, final int height){
	double[][] anchorPoints = null;
	switch (transformation) {
		case 0: {
			anchorPoints = new double[1][3];
			anchorPoints[0][0] = (double)(width / 2);
			anchorPoints[0][1] = (double)(height / 2);
			anchorPoints[0][2] = 1.0;
			break;
		}
		case 1: {
			anchorPoints = new double[3][3];
			anchorPoints[0][0] = (double)(width / 2);
			anchorPoints[0][1] = (double)(height / 2);
			anchorPoints[0][2] = 1.0;
			anchorPoints[1][0] = (double)(width / 2);
			anchorPoints[1][1] = (double)(height / 4);
			anchorPoints[1][2] = 1.0;
			anchorPoints[2][0] = (double)(width / 2);
			anchorPoints[2][1] = (double)((3 * height) / 4);
			anchorPoints[2][2] = 1.0;
			break;
		}
		case 2: {
			anchorPoints = new double[2][3];
			anchorPoints[0][0] = (double)(width / 4);
			anchorPoints[0][1] = (double)(height / 2);
			anchorPoints[0][2] = 1.0;
			anchorPoints[1][0] = (double)((3 * width) / 4);
			anchorPoints[1][1] = (double)(height / 2);
			anchorPoints[1][2] = 1.0;
			break;
		}
		case 3: {
			anchorPoints = new double[3][3];
			anchorPoints[0][0] = (double)(width / 2);
			anchorPoints[0][1] = (double)(height / 4);
			anchorPoints[0][2] = 1.0;
			anchorPoints[1][0] = (double)(width / 4);
			anchorPoints[1][1] = (double)((3 * height) / 4);
			anchorPoints[1][2] = 1.0;
			anchorPoints[2][0] = (double)((3 * width) / 4);
			anchorPoints[2][1] = (double)((3 * height) / 4);
			anchorPoints[2][2] = 1.0;
			break;
		}
		default: {
			IJ.error("Unexpected transformation");
			return null;
		}
	}
	return anchorPoints;
}
/*-------------------------------------------------------------------*/
//return C = A * B
private double[][] matrixMultiply(double[][] A, double[][] B) {
    int mA = A.length;
    int nA = A[0].length;
    int mB = B.length;
    int nB = A[0].length;
    if (nA != mB) throw new RuntimeException("Illegal matrix dimensions.");
    double[][] C = new double[mA][nB];
    for (int i = 0; i < mA; i++)
        for (int j = 0; j < nB; j++)
            for (int k = 0; k < nA; k++)
                C[i][j] += (A[i][k] * B[k][j]);
    return C;
}
/*------------------------------------------------------------------*/
private void computeStatistics (
	final ImageStack imp,
	final double[] average,
	final double[][] scatterMatrix
) {
	int length = imp.getWidth() * imp.getHeight();
	double r;
	double g;
	double b;
	if (imp.getProcessor(1).getPixels() instanceof byte[]) {
		final IndexColorModel icm =
			(IndexColorModel)imp.getProcessor(1).getColorModel();
		final int mapSize = icm.getMapSize();
		final byte[] reds = new byte[mapSize];
		final byte[] greens = new byte[mapSize];
		final byte[] blues = new byte[mapSize];	
		icm.getReds(reds); 
		icm.getGreens(greens); 
		icm.getBlues(blues);
		final double[] histogram = new double[mapSize];
		for (int k = 0; (k < mapSize); k++) {
			histogram[k] = 0.0;
		}
		for (int s = 1; (s <= imp.getSize()); s++) {
			//imp.setSliceWithoutUpdate(s);
			final byte[] pixels = (byte[])imp.getProcessor(s).getPixels();
			for (int k = 0; (k < length); k++) {
				histogram[pixels[k] & 0xFF]++;
			}
		}
		for (int k = 0; (k < mapSize); k++) {
			r = (double)(reds[k] & 0xFF);
			g = (double)(greens[k] & 0xFF);
			b = (double)(blues[k] & 0xFF);
			average[0] += histogram[k] * r;
			average[1] += histogram[k] * g;
			average[2] += histogram[k] * b;
			scatterMatrix[0][0] += histogram[k] * r * r;
			scatterMatrix[0][1] += histogram[k] * r * g;
			scatterMatrix[0][2] += histogram[k] * r * b;
			scatterMatrix[1][1] += histogram[k] * g * g;
			scatterMatrix[1][2] += histogram[k] * g * b;
			scatterMatrix[2][2] += histogram[k] * b * b;
		}
	}
	else if (imp.getProcessor(1).getPixels() instanceof int[]) {
		for (int s = 1; (s <= imp.getSize()); s++) {
			//imp.setSliceWithoutUpdate(s);
			final int[] pixels = (int[])imp.getProcessor(s).getPixels();
			for (int k = 0; (k < length); k++) {
				r = (double)((pixels[k] & 0x00FF0000) >>> 16);
				g = (double)((pixels[k] & 0x0000FF00) >>> 8);
				b = (double)(pixels[k] & 0x000000FF);
				average[0] += r;
				average[1] += g;
				average[2] += b;
				scatterMatrix[0][0] += r * r;
				scatterMatrix[0][1] += r * g;
				scatterMatrix[0][2] += r * b;
				scatterMatrix[1][1] += g * g;
				scatterMatrix[1][2] += g * b;
				scatterMatrix[2][2] += b * b;
			}
		}
	}
	else {
		IJ.error("Internal type mismatch");
	}
	length *= imp.getSize();
	average[0] /= (double)length;
	average[1] /= (double)length;
	average[2] /= (double)length;
	scatterMatrix[0][0] /= (double)length;
	scatterMatrix[0][1] /= (double)length;
	scatterMatrix[0][2] /= (double)length;
	scatterMatrix[1][1] /= (double)length;
	scatterMatrix[1][2] /= (double)length;
	scatterMatrix[2][2] /= (double)length;
	scatterMatrix[0][0] -= average[0] * average[0];
	scatterMatrix[0][1] -= average[0] * average[1];
	scatterMatrix[0][2] -= average[0] * average[2];
	scatterMatrix[1][1] -= average[1] * average[1];
	scatterMatrix[1][2] -= average[1] * average[2];
	scatterMatrix[2][2] -= average[2] * average[2];
	scatterMatrix[2][1] = scatterMatrix[1][2];
	scatterMatrix[2][0] = scatterMatrix[0][2];
	scatterMatrix[1][0] = scatterMatrix[0][1];
} /* computeStatistics */

/*------------------------------------------------------------------*/
private double[] getColorWeightsFromPrincipalComponents (
	final ImageStack imp,
	final int slice
) {
	final double[] average = {0.0, 0.0, 0.0};
	final double[][] scatterMatrix =
		{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	computeStatistics(imp, average, scatterMatrix);
	double[] eigenvalue = getEigenvalues(scatterMatrix);
	if ((eigenvalue[0] * eigenvalue[0] + eigenvalue[1] * eigenvalue[1]
		+ eigenvalue[2] * eigenvalue[2]) <= TINY) {
		return(getLuminanceFromCCIR601());
	}
	double bestEigenvalue = getLargestAbsoluteEigenvalue(eigenvalue);
	double eigenvector[] = getEigenvector(scatterMatrix, bestEigenvalue);
	final double weight = eigenvector[0] + eigenvector[1] + eigenvector[2];
	if (TINY < Math.abs(weight)) {
		eigenvector[0] /= weight;
		eigenvector[1] /= weight;
		eigenvector[2] /= weight;
	}
	return(eigenvector);
} /* getColorWeightsFromPrincipalComponents */

/*------------------------------------------------------------------*/
private double[] getEigenvalues (
	final double[][] scatterMatrix
) {
	final double[] a = {
		scatterMatrix[0][0] * scatterMatrix[1][1] * scatterMatrix[2][2]
			+ 2.0 * scatterMatrix[0][1] * scatterMatrix[1][2]
			* scatterMatrix[2][0]
			- scatterMatrix[0][1] * scatterMatrix[0][1] * scatterMatrix[2][2]
			- scatterMatrix[1][2] * scatterMatrix[1][2] * scatterMatrix[0][0]
			- scatterMatrix[2][0] * scatterMatrix[2][0] * scatterMatrix[1][1],
		scatterMatrix[0][1] * scatterMatrix[0][1]
			+ scatterMatrix[1][2] * scatterMatrix[1][2]
			+ scatterMatrix[2][0] * scatterMatrix[2][0]
			- scatterMatrix[0][0] * scatterMatrix[1][1]
			- scatterMatrix[1][1] * scatterMatrix[2][2]
			- scatterMatrix[2][2] * scatterMatrix[0][0],
		scatterMatrix[0][0] + scatterMatrix[1][1] + scatterMatrix[2][2],
		-1.0
	};
	double[] RealRoot = new double[3];
	double Q = (3.0 * a[1] - a[2] * a[2] / a[3]) / (9.0 * a[3]);
	double R = (a[1] * a[2] - 3.0 * a[0] * a[3]
		- (2.0 / 9.0) * a[2] * a[2] * a[2] / a[3]) / (6.0 * a[3] * a[3]);
	double Det = Q * Q * Q + R * R;
	if (Det < 0.0) {
		Det = 2.0 * Math.sqrt(-Q);
		R /= Math.sqrt(-Q * Q * Q);
		R = (1.0 / 3.0) * Math.acos(R);
		Q = (1.0 / 3.0) * a[2] / a[3];
		RealRoot[0] = Det * Math.cos(R) - Q;
		RealRoot[1] = Det * Math.cos(R + (2.0 / 3.0) * Math.PI) - Q;
		RealRoot[2] = Det * Math.cos(R + (4.0 / 3.0) * Math.PI) - Q;
		if (RealRoot[0] < RealRoot[1]) {
			if (RealRoot[2] < RealRoot[1]) {
				double Swap = RealRoot[1];
				RealRoot[1] = RealRoot[2];
				RealRoot[2] = Swap;
				if (RealRoot[1] < RealRoot[0]) {
					Swap = RealRoot[0];
					RealRoot[0] = RealRoot[1];
					RealRoot[1] = Swap;
				}
			}
		}
		else {
			double Swap = RealRoot[0];
			RealRoot[0] = RealRoot[1];
			RealRoot[1] = Swap;
			if (RealRoot[2] < RealRoot[1]) {
				Swap = RealRoot[1];
				RealRoot[1] = RealRoot[2];
				RealRoot[2] = Swap;
				if (RealRoot[1] < RealRoot[0]) {
					Swap = RealRoot[0];
					RealRoot[0] = RealRoot[1];
					RealRoot[1] = Swap;
				}
			}
		}
	}
	else if (Det == 0.0) {
		final double P = 2.0 * ((R < 0.0) ? (Math.pow(-R, 1.0 / 3.0))
			: (Math.pow(R, 1.0 / 3.0)));
		Q = (1.0 / 3.0) * a[2] / a[3];
		if (P < 0) {
			RealRoot[0] = P - Q;
			RealRoot[1] = -0.5 * P - Q;
			RealRoot[2] = RealRoot[1];
		}
		else {
			RealRoot[0] = -0.5 * P - Q;
			RealRoot[1] = RealRoot[0];
			RealRoot[2] = P - Q;
		}
	}
	else {
		IJ.error("Warning: complex eigenvalue found; ignoring imaginary part.");
		Det = Math.sqrt(Det);
		Q = ((R + Det) < 0.0) ? (-Math.exp((1.0 / 3.0) * Math.log(-R - Det)))
			: (Math.exp((1.0 / 3.0) * Math.log(R + Det)));
		R = Q + ((R < Det) ? (-Math.exp((1.0 / 3.0) * Math.log(Det - R)))
			: (Math.exp((1.0 / 3.0) * Math.log(R - Det))));
		Q = (-1.0 / 3.0) * a[2] / a[3];
		Det = Q + R;
		RealRoot[0] = Q - R / 2.0;
		RealRoot[1] = RealRoot[0];
		RealRoot[2] = RealRoot[1];
		if (Det < RealRoot[0]) {
			RealRoot[0] = Det;
		}
		else {
			RealRoot[2] = Det;
		}
	}
	return(RealRoot);
} /* end getEigenvalues */

/*------------------------------------------------------------------*/
private double[] getEigenvector (
	final double[][] scatterMatrix,
	final double eigenvalue
) {
	final int n = scatterMatrix.length;
	final double[][] matrix = new double[n][n];
	for (int i = 0; (i < n); i++) {
		System.arraycopy(scatterMatrix[i], 0, matrix[i], 0, n);
		matrix[i][i] -= eigenvalue;
	}
	final double[] eigenvector = new double[n];
	double absMax;
	double max;
	double norm;
	for (int i = 0; (i < n); i++) {
		norm = 0.0;
		for (int j = 0; (j < n); j++) {
			norm += matrix[i][j] * matrix[i][j];
		}
		norm = Math.sqrt(norm);
		if (TINY < norm) {
			for (int j = 0; (j < n); j++) {
				matrix[i][j] /= norm;
			}
		}
	}
	for (int j = 0; (j < n); j++) {
		max = matrix[j][j];
		absMax = Math.abs(max);
		int k = j;
		for (int i = j + 1; (i < n); i++) {
			if (absMax < Math.abs(matrix[i][j])) {
				max = matrix[i][j];
				absMax = Math.abs(max);
				k = i;
			}
		}
		if (k != j) {
			final double[] partialLine = new double[n - j];
			System.arraycopy(matrix[j], j, partialLine, 0, n - j);
			System.arraycopy(matrix[k], j, matrix[j], j, n - j);
			System.arraycopy(partialLine, 0, matrix[k], j, n - j);
		}
		if (TINY < absMax) {
			for (k = 0; (k < n); k++) {
				matrix[j][k] /= max;
			}
		}
		for (int i = j + 1; (i < n); i++) {
			max = matrix[i][j];
			for (k = 0; (k < n); k++) {
				matrix[i][k] -= max * matrix[j][k];
			}
		}
	}
	final boolean[] ignore = new boolean[n];
	int valid = n;
	for (int i = 0; (i < n); i++) {
		ignore[i] = false;
		if (Math.abs(matrix[i][i]) < TINY) {
			ignore[i] = true;
			valid--;
			eigenvector[i] = 1.0;
			continue;
		}
		if (TINY < Math.abs(matrix[i][i] - 1.0)) {
			IJ.error("Insufficient accuracy.");
			eigenvector[0] = 0.212671;
			eigenvector[1] = 0.71516;
			eigenvector[2] = 0.072169;
			return(eigenvector);
		}
		norm = 0.0;
		for (int j = 0; (j < i); j++) {
			norm += matrix[i][j] * matrix[i][j];
		}
		for (int j = i + 1; (j < n); j++) {
			norm += matrix[i][j] * matrix[i][j];
		}
		if (Math.sqrt(norm) < TINY) {
			ignore[i] = true;
			valid--;
			eigenvector[i] = 0.0;
			continue;
		}
	}
	if (0 < valid) {
		double[][] reducedMatrix = new double[valid][valid];
		for (int i = 0, u = 0; (i < n); i++) {
			if (!ignore[i]) {
				for (int j = 0, v = 0; (j < n); j++) {
					if (!ignore[j]) {
						reducedMatrix[u][v] = matrix[i][j];
						v++;
					}
				}
				u++;
			}
		}
		double[] reducedEigenvector = new double[valid];
		for (int i = 0, u = 0; (i < n); i++) {
			if (!ignore[i]) {
				for (int j = 0; (j < n); j++) {
					if (ignore[j]) {
						reducedEigenvector[u] -= matrix[i][j] * eigenvector[j];
					}
				}
				u++;
			}
		}
		reducedEigenvector = linearLeastSquares(reducedMatrix,
			reducedEigenvector);
		for (int i = 0, u = 0; (i < n); i++) {
			if (!ignore[i]) {
				eigenvector[i] = reducedEigenvector[u];
				u++;
			}
		}
	}
	norm = 0.0;
	for (int i = 0; (i < n); i++) {
		norm += eigenvector[i] * eigenvector[i];
	}
	norm = Math.sqrt(norm);
	if (Math.sqrt(norm) < TINY) {
		IJ.error("Insufficient accuracy.");
		eigenvector[0] = 0.212671;
		eigenvector[1] = 0.71516;
		eigenvector[2] = 0.072169;
		return(eigenvector);
	}
	absMax = Math.abs(eigenvector[0]);
	valid = 0;
	for (int i = 1; (i < n); i++) {
		max = Math.abs(eigenvector[i]);
		if (absMax < max) {
			absMax = max;
			valid = i;
		}
	}
	norm = (eigenvector[valid] < 0.0) ? (-norm) : (norm);
	for (int i = 0; (i < n); i++) {
		eigenvector[i] /= norm;
	}
	return(eigenvector);
} /* getEigenvector */

/*------------------------------------------------------------------*/
private ImagePlus getGray32 (
	final String title,
	final ImageStack imp,
	final int slice,
	final double[] colorWeights
) {
	final int length = imp.getWidth() * imp.getHeight();
	final ImagePlus gray32 = new ImagePlus(title,
		new FloatProcessor(imp.getWidth(), imp.getHeight()));
	final float[] gray = (float[])gray32.getProcessor().getPixels();
	double r;
	double g;
	double b;
	if (imp.getProcessor(slice).getPixels() instanceof byte[]) {
		final byte[] pixels = (byte[])imp.getProcessor(slice).getPixels();
		final IndexColorModel icm =
			(IndexColorModel)imp.getProcessor(slice).getColorModel();
		final int mapSize = icm.getMapSize();
		final byte[] reds = new byte[mapSize];
		final byte[] greens = new byte[mapSize];
		final byte[] blues = new byte[mapSize];	
		icm.getReds(reds); 
		icm.getGreens(greens); 
		icm.getBlues(blues);
		int index;
		for (int k = 0; (k < length); k++) {
			index = (int)(pixels[k] & 0xFF);
			r = (double)(reds[index] & 0xFF);
			g = (double)(greens[index] & 0xFF);
			b = (double)(blues[index] & 0xFF);
			gray[k] = (float)(colorWeights[0] * r + colorWeights[1] * g
				+ colorWeights[2] * b);
		}
	}
	else if (imp.getProcessor(slice).getPixels() instanceof int[]) {
		final int[] pixels = (int[])imp.getProcessor(slice).getPixels();
		for (int k = 0; (k < length); k++) {
			r = (double)((pixels[k] & 0x00FF0000) >>> 16);
			g = (double)((pixels[k] & 0x0000FF00) >>> 8);
			b = (double)(pixels[k] & 0x000000FF);
			gray[k] = (float)(colorWeights[0] * r + colorWeights[1] * g
				+ colorWeights[2] * b);
		}
	}
	return(gray32);
} /* getGray32 */

/*------------------------------------------------------------------*/
private double getLargestAbsoluteEigenvalue (
	final double[] eigenvalue
) {
	double best = eigenvalue[0];
	for (int k = 1; (k < eigenvalue.length); k++) {
		if (Math.abs(best) < Math.abs(eigenvalue[k])) {
			best = eigenvalue[k];
		}
		if (Math.abs(best) == Math.abs(eigenvalue[k])) {
			if (best < eigenvalue[k]) {
				best = eigenvalue[k];
			}
		}
	}
	return(best);
} /* getLargestAbsoluteEigenvalue */

/*------------------------------------------------------------------*/
private double[] getLuminanceFromCCIR601 (
) {
	double[] weights = {0.299, 0.587, 0.114};
	return(weights);
} /* getLuminanceFromCCIR601 */

/*------------------------------------------------------------------*/
private double[][] getTransformationMatrix (
	final double[][] fromCoord,
	final double[][] toCoord,
	final int transformation
) {
	double[][] matrix = new double[3][3];
	switch (transformation) {
		case 0: {
			matrix[0][0] = 1.0;
			matrix[0][1] = 0.0;
			matrix[0][2] = toCoord[0][0] - fromCoord[0][0];
			matrix[1][0] = 0.0;
			matrix[1][1] = 1.0;
			matrix[1][2] = toCoord[0][1] - fromCoord[0][1];
			break;
		}
		case 1: {
			final double angle = Math.atan2(fromCoord[2][0] - fromCoord[1][0],
				fromCoord[2][1] - fromCoord[1][1])
				- Math.atan2(toCoord[2][0] - toCoord[1][0],
				toCoord[2][1] - toCoord[1][1]);
			final double c = Math.cos(angle);
			final double s = Math.sin(angle);
			matrix[0][0] = c;
			matrix[0][1] = -s;
			matrix[0][2] = toCoord[0][0]
				- c * fromCoord[0][0] + s * fromCoord[0][1];
			matrix[1][0] = s;
			matrix[1][1] = c;
			matrix[1][2] = toCoord[0][1]
				- s * fromCoord[0][0] - c * fromCoord[0][1];
			break;
		}
		case 2: {
			double[][] a = new double[3][3];
			double[] v = new double[3];
			a[0][0] = fromCoord[0][0];
			a[0][1] = fromCoord[0][1];
			a[0][2] = 1.0;
			a[1][0] = fromCoord[1][0];
			a[1][1] = fromCoord[1][1];
			a[1][2] = 1.0;
			a[2][0] = fromCoord[0][1] - fromCoord[1][1] + fromCoord[1][0];
			a[2][1] = fromCoord[1][0] + fromCoord[1][1] - fromCoord[0][0];
			a[2][2] = 1.0;
			invertGauss(a);
			v[0] = toCoord[0][0];
			v[1] = toCoord[1][0];
			v[2] = toCoord[0][1] - toCoord[1][1] + toCoord[1][0];
			for (int i = 0; (i < 3); i++) {
				matrix[0][i] = 0.0;
				for (int j = 0; (j < 3); j++) {
					matrix[0][i] += a[i][j] * v[j];
				}
			}
			v[0] = toCoord[0][1];
			v[1] = toCoord[1][1];
			v[2] = toCoord[1][0] + toCoord[1][1] - toCoord[0][0];
			for (int i = 0; (i < 3); i++) {
				matrix[1][i] = 0.0;
				for (int j = 0; (j < 3); j++) {
					matrix[1][i] += a[i][j] * v[j];
				}
			}
			break;
		}
		case 3: {
			double[][] a = new double[3][3];
			double[] v = new double[3];
			a[0][0] = fromCoord[0][0];
			a[0][1] = fromCoord[0][1];
			a[0][2] = 1.0;
			a[1][0] = fromCoord[1][0];
			a[1][1] = fromCoord[1][1];
			a[1][2] = 1.0;
			a[2][0] = fromCoord[2][0];
			a[2][1] = fromCoord[2][1];
			a[2][2] = 1.0;
			invertGauss(a);
			v[0] = toCoord[0][0];
			v[1] = toCoord[1][0];
			v[2] = toCoord[2][0];
			for (int i = 0; (i < 3); i++) {
				matrix[0][i] = 0.0;
				for (int j = 0; (j < 3); j++) {
					matrix[0][i] += a[i][j] * v[j];
				}
			}
			v[0] = toCoord[0][1];
			v[1] = toCoord[1][1];
			v[2] = toCoord[2][1];
			for (int i = 0; (i < 3); i++) {
				matrix[1][i] = 0.0;
				for (int j = 0; (j < 3); j++) {
					matrix[1][i] += a[i][j] * v[j];
				}
			}
			break;
		}
		default: {
			IJ.error("Unexpected transformation");
		}
	}
	matrix[2][0] = 0.0;
	matrix[2][1] = 0.0;
	matrix[2][2] = 1.0;
	return(matrix);
} /* end getTransformationMatrix */

/*------------------------------------------------------------------*/
private void invertGauss (
	final double[][] matrix
) {
	final int n = matrix.length;
	final double[][] inverse = new double[n][n];
	for (int i = 0; (i < n); i++) {
		double max = matrix[i][0];
		double absMax = Math.abs(max);
		for (int j = 0; (j < n); j++) {
			inverse[i][j] = 0.0;
			if (absMax < Math.abs(matrix[i][j])) {
				max = matrix[i][j];
				absMax = Math.abs(max);
			}
		}
		inverse[i][i] = 1.0 / max;
		for (int j = 0; (j < n); j++) {
			matrix[i][j] /= max;
		}
	}
	for (int j = 0; (j < n); j++) {
		double max = matrix[j][j];
		double absMax = Math.abs(max);
		int k = j;
		for (int i = j + 1; (i < n); i++) {
			if (absMax < Math.abs(matrix[i][j])) {
				max = matrix[i][j];
				absMax = Math.abs(max);
				k = i;
			}
		}
		if (k != j) {
			final double[] partialLine = new double[n - j];
			final double[] fullLine = new double[n];
			System.arraycopy(matrix[j], j, partialLine, 0, n - j);
			System.arraycopy(matrix[k], j, matrix[j], j, n - j);
			System.arraycopy(partialLine, 0, matrix[k], j, n - j);
			System.arraycopy(inverse[j], 0, fullLine, 0, n);
			System.arraycopy(inverse[k], 0, inverse[j], 0, n);
			System.arraycopy(fullLine, 0, inverse[k], 0, n);
		}
		for (k = 0; (k <= j); k++) {
			inverse[j][k] /= max;
		}
		for (k = j + 1; (k < n); k++) {
			matrix[j][k] /= max;
			inverse[j][k] /= max;
		}
		for (int i = j + 1; (i < n); i++) {
			for (k = 0; (k <= j); k++) {
				inverse[i][k] -= matrix[i][j] * inverse[j][k];
			}
			for (k = j + 1; (k < n); k++) {
				matrix[i][k] -= matrix[i][j] * matrix[j][k];
				inverse[i][k] -= matrix[i][j] * inverse[j][k];
			}
		}
	}
	for (int j = n - 1; (1 <= j); j--) {
		for (int i = j - 1; (0 <= i); i--) {
			for (int k = 0; (k <= j); k++) {
				inverse[i][k] -= matrix[i][j] * inverse[j][k];
			}
			for (int k = j + 1; (k < n); k++) {
				matrix[i][k] -= matrix[i][j] * matrix[j][k];
				inverse[i][k] -= matrix[i][j] * inverse[j][k];
			}
		}
	}
	for (int i = 0; (i < n); i++) {
		System.arraycopy(inverse[i], 0, matrix[i], 0, n);
	}
} /* end invertGauss */

/*------------------------------------------------------------------*/
private double[] linearLeastSquares (
	final double[][] A,
	final double[] b
) {
	final int lines = A.length;
	final int columns = A[0].length;
	final double[][] Q = new double[lines][columns];
	final double[][] R = new double[columns][columns];
	final double[] x = new double[columns];
	double s;
	for (int i = 0; (i < lines); i++) {
		for (int j = 0; (j < columns); j++) {
			Q[i][j] = A[i][j];
		}
	}
	QRdecomposition(Q, R);
	for (int i = 0; (i < columns); i++) {
		s = 0.0;
		for (int j = 0; (j < lines); j++) {
			s += Q[j][i] * b[j];
		}
		x[i] = s;
	}
	for (int i = columns - 1; (0 <= i); i--) {
		s = R[i][i];
		if ((s * s) == 0.0) {
			x[i] = 0.0;
		}
		else {
			x[i] /= s;
		}
		for (int j = i - 1; (0 <= j); j--) {
			x[j] -= R[j][i] * x[i];
		}
	}
	return(x);
} /* end linearLeastSquares */

/*------------------------------------------------------------------*/
private void QRdecomposition (
	final double[][] Q,
	final double[][] R
) {
	final int lines = Q.length;
	final int columns = Q[0].length;
	final double[][] A = new double[lines][columns];
	double s;
	for (int j = 0; (j < columns); j++) {
		for (int i = 0; (i < lines); i++) {
			A[i][j] = Q[i][j];
		}
		for (int k = 0; (k < j); k++) {
			s = 0.0;
			for (int i = 0; (i < lines); i++) {
				s += A[i][j] * Q[i][k];
			}
			for (int i = 0; (i < lines); i++) {
				Q[i][j] -= s * Q[i][k];
			}
		}
		s = 0.0;
		for (int i = 0; (i < lines); i++) {
			s += Q[i][j] * Q[i][j];
		}
		if ((s * s) == 0.0) {
			s = 0.0;
		}
		else {
			s = 1.0 / Math.sqrt(s);
		}
		for (int i = 0; (i < lines); i++) {
			Q[i][j] *= s;
		}
	}
	for (int i = 0; (i < columns); i++) {
		for (int j = 0; (j < i); j++) {
			R[i][j] = 0.0;
		}
		for (int j = i; (j < columns); j++) {
			R[i][j] = 0.0;
			for (int k = 0; (k < lines); k++) {
				R[i][j] += Q[k][i] * A[k][j];
			}
		}
	}
} /* end QRdecomposition */

/*------------------------------------------------------------------*/
private void registerSlice (
	final ImagePlus target,
	final ImagePlus imp,
	final ImageStack stack,
	final int width,
	final int height,
	final int transformation,
	final double[][] anchorPoints,
	final double[] colorWeights,
	final int t,
	final int s,
	boolean firstTime
) {
	ImagePlus source=null;
	try {
		Object turboReg = null;
		Object turboRegInv = null;
		Method method = null;
		double[][] sourcePoints = null;
		double[][] targetPoints = null;
		double[][] localTransform = null;
		switch (imp.getType()) {
			case ImagePlus.COLOR_256:
			case ImagePlus.COLOR_RGB: {
				source = getGray32("StackRegSource", stack, s, colorWeights);
				break;
			}
			case ImagePlus.GRAY8: {
				source = new ImagePlus("StackRegSource", new ByteProcessor(
					width, height, (byte[])stack.getProcessor(s).getPixels(),
					stack.getProcessor(s).getColorModel()));
				break;
			}
			case ImagePlus.GRAY16: {
				source = new ImagePlus("StackRegSource", new ShortProcessor(
					width, height, (short[])stack.getProcessor(s).getPixels(),
					stack.getProcessor(s).getColorModel()));
				break;
			}
			case ImagePlus.GRAY32: {
				source = new ImagePlus("StackRegSource", new FloatProcessor(
					width, height, (float[])stack.getProcessor(s).getPixels(),
					stack.getProcessor(s).getColorModel()));
				break;
			}
			default: {
				IJ.error("Unexpected image type");
				return;
			}
		}
		final FileSaver sourceFile = new FileSaver(source);
		final String sourcePathAndFileName =
			IJ.getDirectory("temp") + source.getTitle();
		sourceFile.saveAsTiff(sourcePathAndFileName);
		//IJ.showMessage("MSTStackReg_",sourcePathAndFileName); 
		final FileSaver targetFile = new FileSaver(target);
		final String targetPathAndFileName =
			IJ.getDirectory("temp") + target.getTitle();
		targetFile.saveAsTiff(targetPathAndFileName);
		//IJ.showMessage("MSTStackReg_",targetPathAndFileName); 
		switch (transformation) {
			case 0: {
				turboReg = IJ.runPlugIn("TurboReg_", "-align"
					+ " -file " + sourcePathAndFileName
					+ " 0 0 " + (width - 1) + " " + (height - 1)
					+ " -file " + targetPathAndFileName
					+ " 0 0 " + (width - 1) + " " + (height - 1)
					+ " -translation"
					+ " " + (width / 2) + " " + (height / 2)
					+ " " + (width / 2) + " " + (height / 2)
					+ " -hideOutput"
				);
				if(firstTime)
					turboRegInv = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -translation"
						+ " " + (width / 2) + " " + (height / 2)
						+ " " + (width / 2) + " " + (height / 2)
						+ " -hideOutput"
					);
				break;
			}
			case 1: {
				turboReg = IJ.runPlugIn("TurboReg_", "-align"
					+ " -file " + sourcePathAndFileName
					+ " 0 0 " + (width - 1) + " " + (height - 1)
					+ " -file " + targetPathAndFileName
					+ " 0 0 " + (width - 1) + " " + (height - 1)
					+ " -rigidBody"
					+ " " + (width / 2) + " " + (height / 2)
					+ " " + (width / 2) + " " + (height / 2)
					+ " " + (width / 2) + " " + (height / 4)
					+ " " + (width / 2) + " " + (height / 4)
					+ " " + (width / 2) + " " + ((3 * height) / 4)
					+ " " + (width / 2) + " " + ((3 * height) / 4)
					+ " -hideOutput"
				);
				if(firstTime)
					turboRegInv = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -rigidBody"
						+ " " + (width / 2) + " " + (height / 2)
						+ " " + (width / 2) + " " + (height / 2)
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 2) + " " + ((3 * height) / 4)
						+ " " + (width / 2) + " " + ((3 * height) / 4)
						+ " -hideOutput"
					);
				break;
			}
			case 2: {
				turboReg = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -scaledRotation"
						+ " " + (width / 4) + " " + (height / 2)
						+ " " + (width / 4) + " " + (height / 2)
						+ " " + ((3 * width) / 4) + " " + (height / 2)
						+ " " + ((3 * width) / 4) + " " + (height / 2)
						+ " -hideOutput"
					);
				if(firstTime)
					turboRegInv = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -scaledRotation"
						+ " " + (width / 4) + " " + (height / 2)
						+ " " + (width / 4) + " " + (height / 2)
						+ " " + ((3 * width) / 4) + " " + (height / 2)
						+ " " + ((3 * width) / 4) + " " + (height / 2)
						+ " -hideOutput"
					);
				break;
			}
			case 3: {
				turboReg = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -affine"
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 4) + " " + ((3 * height) / 4)
						+ " " + (width / 4) + " " + ((3 * height) / 4)
						+ " " + ((3 * width) / 4) + " " + ((3 * height) / 4)
						+ " " + ((3 * width) / 4) + " " + ((3 * height) / 4)
						+ " -hideOutput"
					);
				if(firstTime)
					turboRegInv = IJ.runPlugIn("TurboReg_", "-align"
						+ " -file " + targetPathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -file " + sourcePathAndFileName
						+ " 0 0 " + (width - 1) + " " + (height - 1)
						+ " -affine"
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 2) + " " + (height / 4)
						+ " " + (width / 4) + " " + ((3 * height) / 4)
						+ " " + (width / 4) + " " + ((3 * height) / 4)
						+ " " + ((3 * width) / 4) + " " + ((3 * height) / 4)
						+ " " + ((3 * width) / 4) + " " + ((3 * height) / 4)
						+ " -hideOutput"
					);
				break;
			}
			default: {
				IJ.error("Unexpected transformation");
				return;
			}
		}
		if (turboReg == null) {
			throw(new ClassNotFoundException());
		}
		method = turboReg.getClass().getMethod("getSourcePoints",
			(Class[])null);
		sourcePoints = (double[][])method.invoke(turboReg);
		method = turboReg.getClass().getMethod("getTargetPoints",
			(Class[])null);
		targetPoints = (double[][])method.invoke(turboReg);
		localTransform = getTransformationMatrix(targetPoints, sourcePoints,
			transformation);
		method = turboReg.getClass().getMethod("getTransformedImage", (Class[])null);
		ImagePlus transformedImage = (ImagePlus)method.invoke(turboReg);
		double cost = mad(target,transformedImage);
		//IJ.showMessage("MSTStackReg_",IJ.d2s(target.getType(),0)); 
		
		// setting transformation matrix
		for (int i=0; i<3; i++){
			for (int j=0; j<3; j++){
				if(t<s)
					transformationMatrix[t-1][s-t-1][i][j] = localTransform[i][j];
				else
					transformationMatrix[t-1][graphWidth + t-s-1][i][j] = localTransform[i][j];
			}
		}
		// inverse transformation can be computed by registering the target image to the source and is more accurate
		// however rgeitration computation slows down the conmputation. So, following is the faster way to compute inverse transform
		if(firstTime){
		if (turboRegInv == null) {
			throw(new ClassNotFoundException());
		}
		method = turboRegInv.getClass().getMethod("getSourcePoints",
			(Class[])null);
		sourcePoints = (double[][])method.invoke(turboRegInv);
		method = turboRegInv.getClass().getMethod("getTargetPoints",
			(Class[])null);
		targetPoints = (double[][])method.invoke(turboRegInv);
		localTransform = getTransformationMatrix(targetPoints, sourcePoints,
			transformation);
		method = turboRegInv.getClass().getMethod("getTransformedImage", (Class[])null);
		transformedImage = (ImagePlus)method.invoke(turboRegInv);
		double costInv = mad(source,transformedImage);
//		localTransform = getTransformationMatrix(sourcePoints, targetPoints, transformation);
		for (int i=0; i<3; i++){
			for (int j=0; j<3; j++){
				if(t<s)
					transformationMatrix[s-1][graphWidth + s-t-1][i][j] = localTransform[i][j];
				else
					transformationMatrix[s-1][t-s-1][i][j] = localTransform[i][j];
			}
		}
		// create edge and add cost to the graph
		if(t<s){
			DefaultWeightedEdge e1 = graph.addEdge(t, s);
			graph.setEdgeWeight(e1, (cost+costInv)/2); 
		}
		else{
			DefaultWeightedEdge e1 = graph.addEdge(s, t);
			graph.setEdgeWeight(e1, (cost+costInv)/2); 
		}
		}
	} catch (NoSuchMethodException e) {
		IJ.error("Unexpected NoSuchMethodException " + e);
		return;
	} catch (IllegalAccessException e) {
		IJ.error("Unexpected IllegalAccessException " + e);
		return;
	} catch (InvocationTargetException e) {
		IJ.error("Unexpected InvocationTargetException " + e);
		return;
	} catch (ClassNotFoundException e) {
		IJ.error("Please download TurboReg_ from\n"
			+ "http://bigwww.epfl.ch/thevenaz/turboreg/");
		return;
	}			
} /* end registerSlice */

/*------------------------------------------------------------------*/

private ImagePlus transformSlice (
		final ImagePlus imp,
		final int width,
		final int height,
		final int transformation,
		final double[][] globalTransform,
		final double[][] anchorPoints,
		final int s
	) {
		ImagePlus source=null;
		imp.setSlice(s);
		try {
			Object turboReg = null;
			Method method = null;
			double[][] sourcePoints = null;
			double[][] targetPoints = null;
					
			switch (imp.getType()) {
				case ImagePlus.COLOR_256: {
					source = new ImagePlus("StackRegSource", new ByteProcessor(
						width, height, (byte[])imp.getProcessor().getPixels(),
						imp.getProcessor().getColorModel()));
					ImageConverter converter = new ImageConverter(source);
					converter.convertToRGB();
					Object turboRegR = null;
					Object turboRegG = null;
					Object turboRegB = null;
					byte[] r = new byte[width * height];
					byte[] g = new byte[width * height];
					byte[] b = new byte[width * height];
					((ColorProcessor)source.getProcessor()).getRGB(r, g, b);
					final ImagePlus sourceR = new ImagePlus("StackRegSourceR",
						new ByteProcessor(width, height));
					final ImagePlus sourceG = new ImagePlus("StackRegSourceG",
						new ByteProcessor(width, height));
					final ImagePlus sourceB = new ImagePlus("StackRegSourceB",
						new ByteProcessor(width, height));
					sourceR.getProcessor().setPixels(r);
					sourceG.getProcessor().setPixels(g);
					sourceB.getProcessor().setPixels(b);
					ImagePlus transformedSourceR = null;
					ImagePlus transformedSourceG = null;
					ImagePlus transformedSourceB = null;
					final FileSaver sourceFileR = new FileSaver(sourceR);
					final String sourcePathAndFileNameR =
						IJ.getDirectory("temp") + sourceR.getTitle();
					sourceFileR.saveAsTiff(sourcePathAndFileNameR);
					final FileSaver sourceFileG = new FileSaver(sourceG);
					final String sourcePathAndFileNameG =
						IJ.getDirectory("temp") + sourceG.getTitle();
					sourceFileG.saveAsTiff(sourcePathAndFileNameG);
					final FileSaver sourceFileB = new FileSaver(sourceB);
					final String sourcePathAndFileNameB =
						IJ.getDirectory("temp") + sourceB.getTitle();
					sourceFileB.saveAsTiff(sourcePathAndFileNameB);
					switch (transformation) {
						case 0: {
							sourcePoints = new double[1][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 1: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width + " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width + " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						case 2: {
							sourcePoints = new double[2][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 3: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						default: {
							IJ.error("Unexpected transformation");
							return(null);
						}
					}
					method = turboRegR.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceR = (ImagePlus)method.invoke(turboRegR);
					method = turboRegG.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceG = (ImagePlus)method.invoke(turboRegG);
					method = turboRegB.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceB = (ImagePlus)method.invoke(turboRegB);
					transformedSourceR.getStack().deleteLastSlice();
					transformedSourceG.getStack().deleteLastSlice();
					transformedSourceB.getStack().deleteLastSlice();
					transformedSourceR.getProcessor().setMinAndMax(0.0, 255.0);
					transformedSourceG.getProcessor().setMinAndMax(0.0, 255.0);
					transformedSourceB.getProcessor().setMinAndMax(0.0, 255.0);
					ImageConverter converterR =
						new ImageConverter(transformedSourceR);
					ImageConverter converterG =
						new ImageConverter(transformedSourceG);
					ImageConverter converterB =
						new ImageConverter(transformedSourceB);
					converterR.convertToGray8();
					converterG.convertToGray8();
					converterB.convertToGray8();
					final IndexColorModel icm =
						(IndexColorModel)imp.getProcessor().getColorModel();
					final byte[] pixels = (byte[])imp.getProcessor().getPixels();
					r = (byte[])transformedSourceR.getProcessor().getPixels();
					g = (byte[])transformedSourceG.getProcessor().getPixels();
					b = (byte[])transformedSourceB.getProcessor().getPixels();
					final int[] color = new int[4];
					color[3] = 255;
					for (int k = 0; (k < pixels.length); k++) {
						color[0] = (int)(r[k] & 0xFF);
						color[1] = (int)(g[k] & 0xFF);
						color[2] = (int)(b[k] & 0xFF);
						pixels[k] = (byte)icm.getDataElement(color, 0);
					}
					break;
				}
				case ImagePlus.COLOR_RGB: {
					Object turboRegR = null;
					Object turboRegG = null;
					Object turboRegB = null;
					final byte[] r = new byte[width * height];
					final byte[] g = new byte[width * height];
					final byte[] b = new byte[width * height];
					((ColorProcessor)imp.getProcessor()).getRGB(r, g, b);
					final ImagePlus sourceR = new ImagePlus("StackRegSourceR",
						new ByteProcessor(width, height));
					final ImagePlus sourceG = new ImagePlus("StackRegSourceG",
						new ByteProcessor(width, height));
					final ImagePlus sourceB = new ImagePlus("StackRegSourceB",
						new ByteProcessor(width, height));
					sourceR.getProcessor().setPixels(r);
					sourceG.getProcessor().setPixels(g);
					sourceB.getProcessor().setPixels(b);
					ImagePlus transformedSourceR = null;
					ImagePlus transformedSourceG = null;
					ImagePlus transformedSourceB = null;
					final FileSaver sourceFileR = new FileSaver(sourceR);
					final String sourcePathAndFileNameR =
						IJ.getDirectory("temp") + sourceR.getTitle();
					sourceFileR.saveAsTiff(sourcePathAndFileNameR);
					final FileSaver sourceFileG = new FileSaver(sourceG);
					final String sourcePathAndFileNameG =
						IJ.getDirectory("temp") + sourceG.getTitle();
					sourceFileG.saveAsTiff(sourcePathAndFileNameG);
					final FileSaver sourceFileB = new FileSaver(sourceB);
					final String sourcePathAndFileNameB =
						IJ.getDirectory("temp") + sourceB.getTitle();
					sourceFileB.saveAsTiff(sourcePathAndFileNameB);
					switch (transformation) {
						case 0: {
							sourcePoints = new double[1][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 1: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width + " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						case 2: {
							sourcePoints = new double[2][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 3: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboRegR = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameR
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							if (turboRegR == null) {
								throw(new ClassNotFoundException());
							}
							turboRegG = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameG
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							turboRegB = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileNameB
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						default: {
							IJ.error("Unexpected transformation");
							return(null);
						}
					}
					method = turboRegR.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceR = (ImagePlus)method.invoke(turboRegR);
					method = turboRegG.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceG = (ImagePlus)method.invoke(turboRegG);
					method = turboRegB.getClass().getMethod("getTransformedImage",
						(Class[])null);
					transformedSourceB = (ImagePlus)method.invoke(turboRegB);
					transformedSourceR.getStack().deleteLastSlice();
					transformedSourceG.getStack().deleteLastSlice();
					transformedSourceB.getStack().deleteLastSlice();
					transformedSourceR.getProcessor().setMinAndMax(0.0, 255.0);
					transformedSourceG.getProcessor().setMinAndMax(0.0, 255.0);
					transformedSourceB.getProcessor().setMinAndMax(0.0, 255.0);
					ImageConverter converterR =
						new ImageConverter(transformedSourceR);
					ImageConverter converterG =
						new ImageConverter(transformedSourceG);
					ImageConverter converterB =
						new ImageConverter(transformedSourceB);
					converterR.convertToGray8();
					converterG.convertToGray8();
					converterB.convertToGray8();
					((ColorProcessor)imp.getProcessor()).setRGB(
						(byte[])transformedSourceR.getProcessor().getPixels(),
						(byte[])transformedSourceG.getProcessor().getPixels(),
						(byte[])transformedSourceB.getProcessor().getPixels());
					break;
				}
				case ImagePlus.GRAY8:{

					source = new ImagePlus("StackRegSource", new ByteProcessor(
							width, height, (byte[])imp.getProcessor().getPixels(),
							imp.getProcessor().getColorModel()));
					final FileSaver sourceFile = new FileSaver(source);
					final String sourcePathAndFileName =
						IJ.getDirectory("temp") + source.getTitle();
					sourceFile.saveAsTiff(sourcePathAndFileName);

					switch (transformation) {
						case 0: {
							sourcePoints = new double[1][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 1: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						case 2: {
							sourcePoints = new double[2][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 3: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						default: {
							IJ.error("Unexpected transformation");
							return(null);
						}
					}
					if (turboReg == null) {
						throw(new ClassNotFoundException());
					}
					method = turboReg.getClass().getMethod("getTransformedImage",
						(Class[])null);
					ImagePlus transformedSource =
						(ImagePlus)method.invoke(turboReg);
					transformedSource.getStack().deleteLastSlice();
					switch (imp.getType()) {
						case ImagePlus.GRAY8: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 255.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray8();
							break;
						}
						case ImagePlus.GRAY16: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 65535.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray16();
							break;
						}
						case ImagePlus.GRAY32: {
							break;
						}
						default: {
							IJ.error("Unexpected image type");
							return(null);
						}
					}
					imp.setProcessor(null, transformedSource.getProcessor());
					break;
				}
				case ImagePlus.GRAY16:{
					source = new ImagePlus("StackRegSource", new ShortProcessor(
							width, height, (short[])imp.getProcessor().getPixels(),
							imp.getProcessor().getColorModel()));
					final FileSaver sourceFile = new FileSaver(source);
					final String sourcePathAndFileName =
						IJ.getDirectory("temp") + source.getTitle();
					sourceFile.saveAsTiff(sourcePathAndFileName);

					switch (transformation) {
						case 0: {
							sourcePoints = new double[1][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 1: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						case 2: {
							sourcePoints = new double[2][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 3: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						default: {
							IJ.error("Unexpected transformation");
							return(null);
						}
					}
					if (turboReg == null) {
						throw(new ClassNotFoundException());
					}
					method = turboReg.getClass().getMethod("getTransformedImage",
						(Class[])null);
					ImagePlus transformedSource =
						(ImagePlus)method.invoke(turboReg);
					transformedSource.getStack().deleteLastSlice();
					switch (imp.getType()) {
						case ImagePlus.GRAY8: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 255.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray8();
							break;
						}
						case ImagePlus.GRAY16: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 65535.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray16();
							break;
						}
						case ImagePlus.GRAY32: {
							break;
						}
						default: {
							IJ.error("Unexpected image type");
							return(null);
						}
					}
					imp.setProcessor(null, transformedSource.getProcessor());
					break;
				}
				case ImagePlus.GRAY32: {
					source = new ImagePlus("StackRegSource", new FloatProcessor(
							width, height, (float[])imp.getProcessor().getPixels(),
							imp.getProcessor().getColorModel()));
					final FileSaver sourceFile = new FileSaver(source);
					final String sourcePathAndFileName =
						IJ.getDirectory("temp") + source.getTitle();
					sourceFile.saveAsTiff(sourcePathAndFileName);

					switch (transformation) {
						case 0: {
							sourcePoints = new double[1][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -translation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 1: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -rigidBody"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + (width / 2)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						case 2: {
							sourcePoints = new double[2][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -scaledRotation"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 4)
								+ " " + (height / 2)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + ((3 * width) / 4)
								+ " " + (height / 2)
								+ " -hideOutput"
							);
							break;
						}
						case 3: {
							sourcePoints = new double[3][3];
							for (int i = 0; (i < 3); i++) {
								sourcePoints[0][i] = 0.0;
								sourcePoints[1][i] = 0.0;
								sourcePoints[2][i] = 0.0;
								for (int j = 0; (j < 3); j++) {
									sourcePoints[0][i] += globalTransform[i][j]
										* anchorPoints[0][j];
									sourcePoints[1][i] += globalTransform[i][j]
										* anchorPoints[1][j];
									sourcePoints[2][i] += globalTransform[i][j]
										* anchorPoints[2][j];
								}
							}
							turboReg = IJ.runPlugIn("TurboReg_", "-transform"
								+ " -file " + sourcePathAndFileName
								+ " " + width
								+ " " + height
								+ " -affine"
								+ " " + sourcePoints[0][0]
								+ " " + sourcePoints[0][1]
								+ " " + (width / 2)
								+ " " + (height / 4)
								+ " " + sourcePoints[1][0]
								+ " " + sourcePoints[1][1]
								+ " " + (width / 4)
								+ " " + ((3 * height) / 4)
								+ " " + sourcePoints[2][0]
								+ " " + sourcePoints[2][1]
								+ " " + ((3 * width) / 4)
								+ " " + ((3 * height) / 4)
								+ " -hideOutput"
							);
							break;
						}
						default: {
							IJ.error("Unexpected transformation");
							return(null);
						}
					}
					if (turboReg == null) {
						throw(new ClassNotFoundException());
					}
					method = turboReg.getClass().getMethod("getTransformedImage",
						(Class[])null);
					ImagePlus transformedSource =
						(ImagePlus)method.invoke(turboReg);
					transformedSource.getStack().deleteLastSlice();
					switch (imp.getType()) {
						case ImagePlus.GRAY8: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 255.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray8();
							break;
						}
						case ImagePlus.GRAY16: {
							transformedSource.getProcessor().setMinAndMax(
								0.0, 65535.0);
							final ImageConverter converter =
								new ImageConverter(transformedSource);
							converter.convertToGray16();
							break;
						}
						case ImagePlus.GRAY32: {
							break;
						}
						default: {
							IJ.error("Unexpected image type");
							return(null);
						}
					}
					imp.setProcessor(null, transformedSource.getProcessor());
					break;
				}
				default: {
					IJ.error("Unexpected image type");
					return(null);
				}
			}
		} catch (NoSuchMethodException e) {
			IJ.error("Unexpected NoSuchMethodException " + e);
			return(null);
		} catch (IllegalAccessException e) {
			IJ.error("Unexpected IllegalAccessException " + e);
			return(null);
		} catch (InvocationTargetException e) {
			IJ.error("Unexpected InvocationTargetException " + e);
			return(null);
		} catch (ClassNotFoundException e) {
			IJ.error("Please download TurboReg_ from\n"
				+ "http://bigwww.epfl.ch/thevenaz/turboreg/");
			return(null);
		}
		return(source);
	} /* end transformSlice */


} /* end class StackReg_ */