Affine + Consistent Elastic 2D Image Registration

This macro aligns a stack of images taking the first image as reference in two steps: MOPS and bUnwarpJ. The parameters from both methods are introduced in a common dialog.
Requirements

You need to have previously installed the plugins SIFT/MOPS from Stephan Saalfeld and bUnwarpJ from Ignacio Arganda-Carreras.
Input

Stack of 8, 16, 32-bit or RGB-Color images.
Output

Aligned stack.
Masks stack.
Authors

Marta Rivera-Alba
Ignacio Arganda-Carreras
Version and updates

March 11th, 2009: updated to be used with bUnwarpJ 2.5.
November 10th, 2008: changed “divergency_weight” to “divergence_weight” in order to run bUnwarpJ last release.
/**
 * Affine + Consistent Elastic 2D Image Registration macro for ImageJ(C).
 * Copyright (C) 2008 Marta Rivera-Alba and Ignacio Arganda-Carreras 
 *
 * More information at http://biocomp.cnb.csic.es/%7Eiarganda/bUnwarpJ/
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt )
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 * 
 */
 
/**
 * DESCRIPTION: 
 *   This macro aligns a stack of images taking the first image as reference.
 * REQUIREMENTS:
 *   ImageJ with 
 *     bUnwarpJ: http://biocomp.cnb.csic.es/%7Eiarganda/bUnwarpJ/
 *     MOPS: http://pacific.mpi-cbg.de/wiki/index.php/Feature_Extraction
 * INPUT:
 *   Stack of 8, 16, 32-bit or RGB-Color images.
 * OUPUT:
 *   Aligned stack.
 *   Masks stack.
 * AUTHORS:
 *   Marta Rivera-Alba - marta.rivera@uam.es 
 *   Ignacio Arganda-Carreras - ignacio.arganda@uam.es
 * VERSION:
 *   November 10th, 2008.
 *
 */
 
// Create basic dialog
listTransformation = newArray("Translation", "Rigid", "Affine");
listRegistrationMode = newArray("Fast", "Accurate","Mono");
listInitialDef = newArray("Very Coarse","Coarse","Fine","Very Fine");
listFinalDef = newArray("Very Coarse","Coarse","Fine","Very Fine","Super Fine"); 
listImgSub = newArray("0","1","2","3","4","5","6","7");
 
Dialog.create("Affine + Consistent Elastic 2D Registration");
 
Dialog.addMessage("*********** MOPS options ***********");
Dialog.addCheckbox("MOPS landmarks extraction",true)
Dialog.addNumber("steps per scale octave", 3);
Dialog.addNumber("initial gaussian blur", 1.60);
Dialog.addNumber("feature descriptor width", 16);
Dialog.addNumber("minimum image size", 64);
Dialog.addNumber("maximum image size", 1024);
Dialog.addNumber("closest/next closest ratio", 0.92);
Dialog.addNumber("maximal alignment error", 25.0);
Dialog.addNumber("inlier ratio",0.05);
Dialog.addCheckbox("upscale image first",true);
Dialog.addChoice("transformation class",listTransformation, "Affine");
 
Dialog.addMessage("********* bUnwarpJ options *********");
Dialog.addChoice("Registration Mode", listRegistrationMode, "Mono");
Dialog.addChoice("Image_Subsample_Factor", listImgSub, "0");
Dialog.addChoice("Initial Deformation",listInitialDef);
Dialog.addChoice("Final Deformation",listFinalDef, "Fine");
Dialog.addNumber("Divergence Weight", 0.0);
Dialog.addNumber("Curl Weight", 0.0);
Dialog.addNumber("Landmark Weight",0.0);
Dialog.addNumber("Image Weight", 1.0);
Dialog.addNumber("Consistency Weight", 10.0);
Dialog.addNumber("Stop Threshold", 0.01);
Dialog.addCheckbox("Verbose",false);
Dialog.addCheckbox("Save Transformations",false);
 
Dialog.show();
 
 
// Get dialog input values
useMOPS=Dialog.getCheckbox;
spso=Dialog.getNumber;
igb=Dialog.getNumber;
fdw=Dialog.getNumber;
 
minis=Dialog.getNumber;
maxis=Dialog.getNumber;
cncr=Dialog.getNumber;
maxae=Dialog.getNumber;
ir=Dialog.getNumber;
uifbool=Dialog.getCheckbox;
tc=Dialog.getChoice;
 
r=Dialog.getChoice;
img_sub_factor = Dialog.getChoice;
id=Dialog.getChoice;
fd=Dialog.getChoice;
dw=Dialog.getNumber;
curlw=Dialog.getNumber;
lw=Dialog.getNumber;
iw=Dialog.getNumber;
conw=Dialog.getNumber;
st=Dialog.getNumber;
verbosebool=Dialog.getCheckbox;
strbool=Dialog.getCheckbox;
 
// Transform values into macro parameters
if (uifbool==1){uif="upscale_image_first";} else {uif="";}
 
if (id=="Very Coarse")
	{id="[Very Coarse]";} 
else if (id=="Very Fine")
	{id="[Very Fine]";}
 
 
if (fd=="Very Coarse")
	{fd="[Very Coarse]";} 
else if (fd=="Very Fine")
	{fd="[Very Fine]";}
else if (fd=="Super Fine")
	{fd="[Super Fine]";}
 
 
if (verbosebool==1){verbose="verbose";} else {verbose="";}
 
if (strbool==1){str="save_transformations";} else {str="";}
 
 
 
// Get input stack name
inputStackName = getTitle();
 
// Select input stack
selectWindow(inputStackName);
numOfInputSlices = nSlices;
 
// Duplicate stack to store results
run("Duplicate...", "title=[Registration results] duplicate");
// Grayscale images are converted to 32-bit.
if(bitDepth() != 24)
    run("32-bit");
// Duplicate stack to store masks
run("Duplicate...", "title=[Registration masks] duplicate");
// Grayscale images are converted to 32-bit.
if(bitDepth() != 24)
    run("32-bit");
// The two initial masks are all white
setForegroundColor(255, 255, 255);
for(i = 1; i < 3; i++)
{
    setSlice(i);
    run("Select All");
    run("Fill", "slice");
}
 
// Unidirectional registration ("Mono" mode)
if(r == "Mono")
{
    for(i = 1; i < numOfInputSlices; i++)
    {
        sourceName = "source" + i;
        targetName = "target" + i;
        // Take 2 consecutive slices (i and i+1).
        selectWindow("Registration results");
        setSlice(i);
        run("Duplicate...", "title=" + sourceName);    
        selectWindow("Registration results");
        setSlice(i+1);
        run("Duplicate...", "title=" + targetName);
        // Add corresponding mask.
        selectWindow("Registration masks");
        setSlice(i);
        run("Copy");
        selectWindow(sourceName);
        run("Add Slice");
        run("Paste");
        setSlice(1);
 
        if (useMOPS==1)
        {   
 
            // Optional previous MOPS landmarks extraction
            run("Extract MOPS Correspondences", "source_image=" + sourceName + " target_image=" + targetName + " steps_per_scale_octave=" + spso + " initial_gaussian_blur=" + igb + " feature_descriptor_width=" + fdw + " minimum_image_size=" + minis + " maximum_image_size=" + maxis + " closest/next_closest_ratio=" + cncr + " maximal_alignment_error=" + maxae + " inlier_ratio=" + ir + " " +  uif + " transformation_class="+ tc);
 
            print("Finished MOPS between slice " + i + " and " + (i+1));
        }
 
        // Register with bUnwarpJ
        run("bUnwarpJ", "source_image=" + targetName + " target_image=" + sourceName + " registration=" + r + " image_subsample_factor="+ img_sub_factor + " initial_deformation=" + id + " final_deformation=" + fd + " divergence_weight=" + dw + " curl_weight=" + curlw + " landmark_weight="+ lw + " image_weight=" + iw + " consistency_weight=" + conw + " stop_threshold=" + st + " " + verbose + str);        
 
        // Wait bUnwarpJ to finish
        while (isOpen("Registered Source Image") != 1)
        { 
            wait(1000);
        }     
        wait(1500);
        print("Finished MONO bUnwarpJ between slice " + i + " (" + sourceName + ") and " + (i+1) + " (" + targetName + ")");
 
        // Copy registration result ...
        selectWindow("Registered Source Image");
        run("Copy");
 
        // ... into the registration results stack 
        selectWindow("Registration results");
        run("Paste");
 
        // And copy mask too.
        selectWindow("Registered Source Image");
        setSlice(3);
        run("Copy");
        close();
        selectWindow("Registration masks");
        setSlice(i+1);
        run("Paste");
 
        // Close temporary source and target windows
        selectWindow(sourceName);
        close();
        selectWindow(targetName);
        close();
    }
}
else // Bidirectional registration ("Fast" or "Accurate" modes)
{
    for(i = 1; i < numOfInputSlices; i++)
    {
        sourceName = "source" + i;
        targetName = "target" + i;
        // Take 2 consecutive slices (i and i+1).
        selectWindow("Registration results");
        setSlice(i);
        run("Duplicate...", "title=" + sourceName);    
        selectWindow("Registration results");
        setSlice(i+1);
        run("Duplicate...", "title=" + targetName);
        // Add corresponding mask.
        selectWindow("Registration masks");
        setSlice(i);
        run("Copy");
        selectWindow(sourceName);
        run("Add Slice");
        run("Paste");
        setSlice(1);
 
        if (useMOPS==1)
        {   
 
            // Optional previous MOPS landmarks extraction
            run("Extract MOPS Correspondences", "source_image=" + sourceName + " target_image=" + targetName + " steps_per_scale_octave=" + spso + " initial_gaussian_blur=" + igb + " feature_descriptor_width=" + fdw + " minimum_image_size=" + minis + " maximum_image_size=" + maxis + " closest/next_closest_ratio=" + cncr + " maximal_alignment_error=" + maxae + " inlier_ratio=" + ir + " " +  uif + " transformation_class="+ tc);
 
            print("Finished MOPS between slice " + i + " and " + (i+1));
        }
 
        // Register with bUnwarpJ
        run("bUnwarpJ", "source_image=" + sourceName + " target_image=" + targetName + " registration=" + r + " initial_deformation=" + id + " final_deformation=" + fd + " divergence_weight=" + dw + " curl_weight=" + curlw + " landmark_weight="+ lw + " image_weight=" + iw + " consistency_weight=" + conw + " stop_threshold=" + st + " " + verbose + str);        
 
        // Wait bUnwarpJ to finish
        while (isOpen("Registered Source Image") != 1)
        { 
            wait(1000);
        }     
        wait(1500);
        print("Finished bUnwarpJ between slice " + i + " and " + (i+1));
 
        // Close direct registration results window
        selectWindow("Registered Source Image");
        close();
 
        // Copy inverse registration result ...
        selectWindow("Registered Target Image");
        run("Copy");
 
        // ... into the registration results stack 
        selectWindow("Registration results");
        run("Paste");
 
        // And copy mask too.
        selectWindow("Registered Target Image");
        setSlice(3);
        run("Copy");
        close();
        selectWindow("Registration masks");
        setSlice(i+1);
        run("Paste");
 
        // Close temporary source and target windows
        selectWindow(sourceName);
        close();
        selectWindow(targetName);
        close();
    }
}