/*******************************************************************************
 * This macro example performs the registration of two images that are related
 * by a rotation. It checks that the registration is successful by computing the
 * residual difference between the target image and the registered image.
 * i) Prepare the target image.
 * ii) Prepare the source image.
 * iii) Perform the registration.
 * iv) Compute the residual differences and compare to ImageJ.
 ***/

/*******************************************************************************
 * i) Preparation of the target image.
 *
 * 1) Open one of the standard images.
 * 2) Change its type. The reason for this step will become apparent later.
 * 3) Save it in a temporary folder.
 ***/

// 1) Open the Blobs image and set an appropriate name.
run("Blobs (25K)");
unrotatedBlobsTitle = "unrotatedBlobs";
rename(unrotatedBlobsTitle);

// 2) Update the lookup table so that low numeric values are dark and so that
// high numeric values are light, and convert the image to 32-bit. Normalize
// the displayed intensities to the original range.
run("Grays");
run("32-bit");
setMinAndMax(0, 255);

// 3) Save the target image into the temporary directory of ImageJ. This
// illustrates how to deal with a file. We shall illustrate with the source
// image how to refer to an image by the name of its window.
path = getDirectory("temp");
unrotatedBlobsTitle = unrotatedBlobsTitle + ".tif";
run("Tiff...", "save=" + path + unrotatedBlobsTitle);

/*******************************************************************************
 * ii) Preparation of the source image.
 *
 * 1) Open the same standard image and change its type.
 * 2) Perform a 10-degree rotation.
 * 3) Adjoin a mask [this step is optional].
 ***/

// 1) Create another image with the appropriate type and name.
rotatedBlobsTitle = "rotatedBlobs";
selectWindow(unrotatedBlobsTitle);
run("Duplicate...", "title=" + rotatedBlobsTitle);

// 2) The transformation we apply is standard in ImageJ. It is a simple rotation
// by 10 degree.
run("Arbitrarily...", "angle=10 interpolate");

// 3) The next few steps illustrate one way to build a mask. They are optional.
// Their goal is to mask out the invalid parts of the image that appeared in
// reason of the rotation.
width = getWidth();
height = getHeight();
run("Add Slice"); // Create the slice that is holding the mask.
run("Select All"); // Initially, nothing is masked out.
// Filled areas are set to 255 and deleted areas to 0.
run("Colors...", "foreground=white background=black selection=yellow");
run("Fill", "slice");
run("Set Slice...", "slice=1"); // Select the slice that is holding the data.
setTool(8); // We are now going to use the magic wand.
doWand(0, 0); // Pick the top-left corner.
run("Next Slice [>]"); // Transfer the selection to the slice holding the mask.
run("Clear", "slice"); // Set to 0 the area that needs to be masked out.
run("Previous Slice [<]"); // Repeat for the other three corners.
doWand(0, height - 1); // Bottom-left.
run("Next Slice [>]"); run("Clear", "slice"); run("Previous Slice [<]");
doWand(width - 1, height - 1); // Bottom-right.
run("Next Slice [>]"); run("Clear", "slice"); run("Previous Slice [<]");
doWand(width - 1, 0); // Top-right.
run("Next Slice [>]"); run("Clear", "slice"); run("Previous Slice [<]");
run("Select None"); // Get rid of the ROI.
setMinAndMax(0, 255); // Normalize the displayed intensities.

/*******************************************************************************
 * iii) Registration.
 *
 * 1) Call TurboReg.
 * 2) Give an appropriate name and normalize the displayed intensities.
 * 3) Compute and display the rotation angle as recovered by TurboReg.
 ***/

// 1) Call the TurboReg plugin. This is the central element of this example.
run("TurboReg ",
	"-align " // Register the two images that we have just prepared.
	+ "-window " + rotatedBlobsTitle + " "// Source (window reference).
	+ "0 0 " + (width - 1) + " " + (height - 1) + " " // No cropping.
	+ "-file " + path + unrotatedBlobsTitle + " "// Target (file reference).
	+ "0 0 " + (width - 1) + " " + (height - 1) + " " // No cropping.
	+ "-rigidBody " // This corresponds to rotation and translation.
	+ (width / 2) + " " + (height / 2) + " " // Source translation landmark.
	+ (width / 2) + " " + (height / 2) + " " // Target translation landmark.
	+ "0 " + (height / 2) + " " // Source first rotation landmark.
	+ "0 " + (height / 2) + " " // Target first rotation landmark.
	+ (width - 1) + " " + (height / 2) + " " // Source second rotation landmark.
	+ (width - 1) + " " + (height / 2) + " " // Target second rotation landmark.
	+ "-showOutput"); // In case -hideOutput is selected, the only way to
	// retrieve the registration results is by the way of another plugin.

// 2) Give a more memorable name, and a better visual aspect.
selectWindow("Output");
rename("registered");
setMinAndMax(0, 255);

// 3) The numeric results are stored inside ImageJ's results table. The first
// line (index [0]) corresponds to the coordinates of the pair of landmarks that
// determine the translation component. This translation might be non-zero
// because of geometric dicrepancies between ImageJ's coordinate system and
// TurboReg's. Together, the two remaining lines determine the rotation: we can
// compute the rotation angle by using simple trigonometric relations. Because
// we are undoing the rotation of step ii), the sign of the rotation angle is
// opposite to that we used in step ii).
sourceX0 = getResult("sourceX", 0); // First line of the table.
sourceY0 = getResult("sourceY", 0);
targetX0 = getResult("targetX", 0);
targetY0 = getResult("targetY", 0);
sourceX1 = getResult("sourceX", 1); // Second line of the table.
sourceY1 = getResult("sourceY", 1);
targetX1 = getResult("targetX", 1);
targetY1 = getResult("targetY", 1);
sourceX2 = getResult("sourceX", 2); // Third line of the table.
sourceY2 = getResult("sourceY", 2);
targetX2 = getResult("targetX", 2);
targetY2 = getResult("targetY", 2);
dx = targetX0 - sourceX0;
dy = targetY0 - sourceY0;
translation = sqrt(dx * dx + dy * dy); // Amount of translation, in pixel units.
dx = sourceX2 - sourceX1;
dy = sourceY2 - sourceY1;
sourceAngle = atan2(dy, dx);
dx = targetX2 - targetX1;
dy = targetY2 - targetY1;
targetAngle = atan2(dy, dx);
rotation = targetAngle - sourceAngle; // Amount of rotation, in radian units.
print("Amount of translation [pixel]: " + translation);
print("Amount of rotation [degree]: " + (rotation * 180.0 / PI));

/*******************************************************************************
 * iv) Residual difference.
 * Let B be the blobs image. Then, the strategy is to compute:
 * C = rotate(B, +10 degree) // Do some rotation.
 * D = register(target = B, source = C) // Try to blindly undo the rotation.
 * E = B - D // Determine the residual registration error.
 * F = rotate(C, -10 degree) // We should have F = B.
 * G = B - F // Residual error due to interpolation alone, not to registration.
 *
 * Note: C and D have already been computed in ii) and iii), respectively.
 * 1) Compute E.
 * 2) Compute F.
 * 3) Compute G.
 * 4) Present all images.
 * 5) Analyze and compare E to G.
 * 6) Conclude.
 ***/

// 1) The main reason for changing the image type is to allow the computation of
// a difference with the image resulting from TurboReg, which is always of type
// 32-bit. ImageJ enforces that differences can be computed on images of same
// type only.
run("Image Calculator...",
	"image1=" + unrotatedBlobsTitle + " "
	+ "operation=Subtract " // Compute the difference E = B - D.
	+ "image2=registered create");
selectWindow("Result of " + unrotatedBlobsTitle);
rename("diffAfterRegistration"); // Give a descriptive name.
setMinAndMax(-8, 8); // Visual formatting.

// 2) Undo the rotation of C to get F = backAndForth.
selectWindow(rotatedBlobsTitle);
run("Duplicate...", "title=backAndForth");
run("Arbitrarily...", "angle=-10 interpolate"); // Undo the rotation.

// 3) If the rotation would be perfect, the images B (unrotatedBlobs) and F
// (backAndForth) would be identical. This is not the case because of
// interpolation artefacts. We check now by how much F differs from B.
run("Image Calculator...",
	"image1=" + unrotatedBlobsTitle + " "
	+ "operation=Subtract " // Compute the difference G = B - F.
	+ "image2=backAndForth create");
selectWindow("Result of " + unrotatedBlobsTitle);
rename("diffAfterBackAndForth"); // Give a descriptive name.
setMinAndMax(-8, 8); // Visual formatting.

// 4) Present all images
run("Tile");

// 5) To check that the registration did succeed, we compute the standard
// deviation of the residual error.
selectWindow("diffAfterRegistration");
makeOval(2, 2, width - 4, height - 4); // Look away from the boundaries.
run("Histogram", "number=400 histogram_min=-20 histogram_max=20");
selectWindow("diffAfterBackAndForth");
makeOval(2, 2, width - 4, height - 4); // Look away from the boundaries.
run("Histogram", "number=400 histogram_min=-20 histogram_max=20");

// 6) The visual inspection of diffAfterRegistration reveals that the
// registration has been successful--often, a mildly successful registration
// would reveal itself by ``shadows'' in the difference image (which is not the
// case here), while a failed registration results in high values of the
// standard deviation (which is not the case either). Moreover, since the
// standard deviation of the difference is lower in the TurboReg case than in
// the back-and-forth case, TurboReg does even better than the best that ImageJ
// can offer! Inspection of the minimal and maximal error confirms that
// statement, which can be explained by a better interpolation model in TurboReg
// (cubic splines) than is available from ImageJ.