package edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration;
/*====================================================================
| Version: July 8, 2006 by Steve Bryson and Carlos �scar
| modified to replace Tiff image output with high-quality Jpeg output
\===================================================================*/
/*====================================================================
| Carlos Oscar Sanchez Sorzano
| Unidad de Bioinformatica
| Centro Nacional de Biotecnologia (CSIC)
| Campus Universidad Autonoma (Cantoblanco)
| E-28049 Madrid
| Spain
|
| phone (CET): +34(91)585.45.10
| fax: +34(91)585.45.06
| RFC-822: coss@cnb.uam.es
| URL: http://biocomp.cnb.uam.es/
\===================================================================*/
/*====================================================================
| This work is based on the following paper:
|
| C.O.S. Sorzano, P. Thevenaz, M. Unser
| Elastic Registration of Biological Images Using Vector-Spline Regularization
| IEEE Transactions on Biomedical Imaging, 52: 652-663 (2005)
|
| This paper is available on-line at
| http://bigwww.epfl.ch/publications/sorzano0501.html
|
| Other relevant on-line publications are available at
| http://bigwww.epfl.ch/publications/
\===================================================================*/
/*====================================================================
| Additional help available at http://bigwww.epfl.ch/thevenaz/UnwarpJ/
|
| You'll be free to use this software for research purposes, but you
| should not redistribute it without our consent. In addition, we expect
| you to include a citation or acknowledgment whenever you present or
| publish results that are based on it.
\===================================================================*/
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.Macro;
import ij.WindowManager;
import ij.gui.GUI;
import ij.gui.ImageCanvas;
import ij.gui.ImageWindow;
import ij.gui.Roi;
import ij.gui.Toolbar;
import ij.io.FileSaver;
import ij.io.Opener;
import ij.measure.Calibration;
import ij.plugin.JpegWriter;
import ij.plugin.PlugIn;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageConverter;
import ij.process.ImageProcessor;
import ij.process.ShortProcessor;
import java.awt.BorderLayout;
import java.awt.Button;
import java.awt.Canvas;
import java.awt.Checkbox;
import java.awt.CheckboxGroup;
import java.awt.Choice;
import java.awt.Color;
import java.awt.Component;
import java.awt.Container;
import java.awt.Dialog;
import java.awt.Event;
import java.awt.FileDialog;
import java.awt.FlowLayout;
import java.awt.Font;
import java.awt.Frame;
import java.awt.Graphics;
import java.awt.GridLayout;
import java.awt.Insets;
import java.awt.Label;
import java.awt.Panel;
import java.awt.Point;
import java.awt.Polygon;
import java.awt.Rectangle;
import java.awt.TextArea;
import java.awt.TextField;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.awt.event.ItemEvent;
import java.awt.event.ItemListener;
import java.awt.event.KeyEvent;
import java.awt.event.KeyListener;
import java.awt.event.MouseEvent;
import java.awt.event.MouseListener;
import java.awt.event.MouseMotionListener;
import java.awt.event.TextEvent;
import java.awt.event.TextListener;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;
import java.util.Stack;
import java.util.StringTokenizer;
import java.util.Vector;
/*====================================================================
| UnwarpJ_
\===================================================================*/
/* Main class.
This is the one called by ImageJ */
/*------------------------------------------------------------------*/
public class UnwarpJ_
implements
PlugIn
{ /* begin class UnwarpJ_ */
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void run (
final String commandLine
) {
String options = Macro.getOptions();
if (!commandLine.equals("")) options = commandLine;
if (options == null) {
Runtime.getRuntime().gc();
final ImagePlus[] imageList = createImageList();
if (imageList.length < 2) {
IJ.error("At least two images are required (stack of color images disallowed)");
return;
}
final unwarpJDialog dialog = new unwarpJDialog(IJ.getInstance(), imageList);
GUI.center(dialog);
dialog.setVisible(true);
} else {
final String[] args = getTokens(options);
if (args.length<1) {
dumpSyntax();
return;
} else {
if (args[0].equals("-help")) dumpSyntax();
else if (args[0].equals("-align")) alignImagesMacro(args);
else if (args[0].equals("-transform")) transformImageMacro(args);
}
return;
}
} /* end run */
/*------------------------------------------------------------------*/
public static void main(String args[]) {
if (args.length<1) {
dumpSyntax();
System.exit(1);
} else {
if (args[0].equals("-help")) dumpSyntax();
else if (args[0].equals("-align")) alignImagesMacro(args);
else if (args[0].equals("-transform")) transformImageMacro(args);
}
System.exit(0);
}
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private static void alignImagesMacro(String args[]) {
if(args.length < 11)
{
dumpSyntax();
System.exit(0);
}
// Check if -output_jpg at the end
int args_len=args.length;
String last_argument=args[args_len-1];
boolean jpeg_output=last_argument.equals("-jpg_output");
if (jpeg_output) args_len--;
// Read input parameters
String fn_target=args[1];
String fn_target_mask=args[2];
String fn_source=args[3];
String fn_source_mask=args[4];
int min_scale_deformation=((Integer) new Integer(args[5])).intValue();
int max_scale_deformation=((Integer) new Integer(args[6])).intValue();
double divWeight=((Double) new Double(args[7])).doubleValue();
double curlWeight=((Double) new Double(args[8])).doubleValue();
double imageWeight=((Double) new Double(args[9])).doubleValue();
String fn_out=args[10];
double landmarkWeight=0;
String fn_landmark="";
if (args_len>=13) {
landmarkWeight=((Double) new Double(args[11])).doubleValue();
fn_landmark=args[12];
}
double stopThreshold=1e-2;
if (args_len>=14)
stopThreshold=((Double) new Double(args[13])).doubleValue();
// Show parameters
IJ.write("Target image : "+fn_target);
IJ.write("Target mask : "+fn_target_mask);
IJ.write("Source image : "+fn_source);
IJ.write("Source mask : "+fn_source_mask);
IJ.write("Min. Scale Deformation : "+min_scale_deformation);
IJ.write("Max. Scale Deformation : "+max_scale_deformation);
IJ.write("Div. Weight : "+divWeight);
IJ.write("Curl Weight : "+curlWeight);
IJ.write("Image Weight : "+imageWeight);
IJ.write("Output: : "+fn_out);
IJ.write("Landmark Weight : "+landmarkWeight);
IJ.write("Landmark file : "+fn_landmark);
IJ.write("JPEG Output : "+jpeg_output);
// Produce side information
int imagePyramidDepth=max_scale_deformation-min_scale_deformation+1;
int min_scale_image=0;
int outputLevel=-1;
//int outputLevel= 1;
boolean showMarquardtOptim=false;
int accurate_mode=1;
boolean saveTransf=true;
String fn_tnf="";
int dot = fn_out.lastIndexOf('.');
if (dot == -1) fn_tnf=fn_out + "_transf.txt";
else fn_tnf=fn_out.substring(0, dot)+"_transf.txt";
// Open target
Opener opener=new Opener();
ImagePlus targetImp;
targetImp=opener.openImage(fn_target);
unwarpJImageModel target =
new unwarpJImageModel(targetImp.getProcessor(), true);
target.setPyramidDepth(imagePyramidDepth+min_scale_image);
target.getThread().start();
unwarpJMask targetMsk =
new unwarpJMask(targetImp.getProcessor(),false);
if (fn_target_mask.equalsIgnoreCase(new String("NULL")) == false)
targetMsk.readFile(fn_target_mask);
unwarpJPointHandler targetPh=null;
// Open source
ImagePlus sourceImp;
sourceImp=opener.openImage(fn_source);
unwarpJImageModel source =
new unwarpJImageModel(sourceImp.getProcessor(), false);
source.setPyramidDepth(imagePyramidDepth+min_scale_image);
source.getThread().start();
unwarpJMask sourceMsk =
new unwarpJMask(sourceImp.getProcessor(),false);
if (fn_source_mask.equalsIgnoreCase(new String("NULL")) == false)
sourceMsk.readFile(fn_source_mask);
unwarpJPointHandler sourcePh=null;
// Load landmarks
if (fn_landmark!="") {
Stack sourceStack = new Stack();
Stack targetStack = new Stack();
unwarpJMiscTools.loadPoints(fn_landmark,sourceStack,targetStack);
sourcePh = new unwarpJPointHandler(sourceImp);
targetPh = new unwarpJPointHandler(targetImp);
while ((!sourceStack.empty()) && (!targetStack.empty())) {
Point sourcePoint = (Point)sourceStack.pop();
Point targetPoint = (Point)targetStack.pop();
sourcePh.addPoint(sourcePoint.x, sourcePoint.y);
targetPh.addPoint(targetPoint.x, targetPoint.y);
}
}
// Join threads
try {
source.getThread().join();
target.getThread().join();
} catch (InterruptedException e) {
IJ.error("Unexpected interruption exception " + e);
}
// Perform registration
ImagePlus output_ip=new ImagePlus();
unwarpJDialog dialog=null;
final unwarpJTransformation warp = new unwarpJTransformation(
sourceImp, targetImp, source, target, sourcePh, targetPh,
sourceMsk, targetMsk, min_scale_deformation, max_scale_deformation,
min_scale_image, divWeight, curlWeight, landmarkWeight, imageWeight,
stopThreshold, outputLevel, showMarquardtOptim, accurate_mode,
saveTransf, fn_tnf, output_ip, dialog);
warp.doRegistration();
// Save results
ImageConverter converter=new ImageConverter(output_ip);
converter.convertToGray16();
FileSaver fs=new FileSaver(output_ip);
if (jpeg_output) {
JpegWriter js = new JpegWriter();
js.setQuality(100);
WindowManager.setTempCurrentImage(output_ip);
js.run(fn_out);
} else
fs.saveAsTiff(fn_out);
}
/*------------------------------------------------------------------*/
private ImagePlus[] createImageList (
) {
final int[] windowList = WindowManager.getIDList();
final Stack stack = new Stack();
for (int k = 0; ((windowList != null) && (k < windowList.length)); k++) {
final ImagePlus imp = WindowManager.getImage(windowList[k]);
final int inputType = imp.getType();
if ((imp.getStackSize() == 1) || (inputType == imp.GRAY8) || (inputType == imp.GRAY16)
|| (inputType == imp.GRAY32)) {
stack.push(imp);
}
}
final ImagePlus[] imageList = new ImagePlus[stack.size()];
int k = 0;
while (!stack.isEmpty()) {
imageList[k++] = (ImagePlus)stack.pop();
}
return(imageList);
} /* end createImageList */
/*------------------------------------------------------------------*/
private static void dumpSyntax (
) {
IJ.write("Purpose: Elastic registration of two images.");
IJ.write(" ");
IJ.write("Usage: unwarpj ");
IJ.write(" -help : SHOWS THIS MESSAGE");
IJ.write("");
IJ.write(" -align : ALIGN TWO IMAGES");
IJ.write(" target_image : In any image format");
IJ.write(" target_mask : In any image format");
IJ.write(" source_image : In any image format");
IJ.write(" source_mask : In any image format");
IJ.write(" min_scale_def : Scale of the coarsest deformation");
IJ.write(" 0 is the coarsest possible");
IJ.write(" max_scale_def : Scale of the finest deformation");
IJ.write(" Div_weight : Weight of the divergence term");
IJ.write(" Curl_weight : Weight of the curl term");
IJ.write(" Image_weight : Weight of the image term");
IJ.write(" Output image : Output result in JPG (100%)");
IJ.write(" Optional parameters :");
IJ.write(" Landmark_weight : Weight of the landmarks");
IJ.write(" Landmark_file : Landmark file");
IJ.write(" stopThreshold : By default 1e-2");
IJ.write("");
IJ.write(" -transform : TRANSFORM AN IMAGE WITH A GIVEN DEFORMATION");
IJ.write(" target_image : In any image format");
IJ.write(" source_image : In any image format");
IJ.write(" transformation_file : As saved by UnwarpJ");
IJ.write(" Output image : Output result in JPG (100%)");
IJ.write("");
IJ.write("Examples:");
IJ.write("Align two images without landmarks and without mask");
IJ.write(" unwarpj -align target.jpg NULL source.jpg NULL 0 2 1 1 1 output.tif");
IJ.write("Align two images with landmarks and mask");
IJ.write(" unwarpj -align target.tif target_mask.tif source.tif source_mask.tif 0 2 1 1 1 1 landmarks.txt output.tif");
IJ.write("Align two images using only landmarks");
IJ.write(" unwarpj -align target.jpg NULL source.jpg NULL 0 2 1 1 0 output.tif 1 landmarks.txt");
IJ.write("Transform the source image with a previously computed transformation");
IJ.write(" unwarpj -transform target.jpg source.jpg transformation.txt output.tif");
IJ.write("");
IJ.write("JPEG Output:");
IJ.write(" If you want to produce JPEG output simply add -jpg_output as the last argument");
IJ.write(" of the alignment or transformation command. For instance:");
IJ.write(" unwarpj -align target.jpg NULL source.jpg NULL 0 2 1 1 1 output.jpg -jpg_output");
IJ.write(" unwarpj -align target.tif target_mask.tif source.tif source_mask.tif 0 2 1 1 1 1 landmarks.txt output.jpg -jpg_output");
IJ.write(" unwarpj -align target.jpg NULL source.jpg NULL 0 2 1 1 0 output.jpg 1 landmarks.txt -jpg_output");
IJ.write(" unwarpj -transform target.jpg source.jpg transformation.txt output.jpg -jpg_output");
} /* end dumpSyntax */
/*------------------------------------------------------------------*/
private String[] getTokens (
final String options
) {
StringTokenizer t = new StringTokenizer(options);
String[] token = new String[t.countTokens()];
for (int k = 0; (k < token.length); k++) {
token[k] = t.nextToken();
}
return(token);
} /* end getTokens */
/*------------------------------------------------------------------*/
private static void transformImageMacro(String args[]) {
// Read input parameters
String fn_target=args[1];
String fn_source=args[2];
String fn_tnf =args[3];
String fn_out =args[4];
// Jpeg output
String last_argument=args[args.length-1];
boolean jpeg_output=last_argument.equals("-jpg_output");
// Show parameters
IJ.write("Target image : "+fn_target);
IJ.write("Source image : "+fn_source);
IJ.write("Transformation file : "+fn_tnf);
IJ.write("Output : "+fn_out);
IJ.write("JPEG output : "+jpeg_output);
// Open target
Opener opener=new Opener();
ImagePlus targetImp;
targetImp=opener.openImage(fn_target);
// Open source
ImagePlus sourceImp;
sourceImp=opener.openImage(fn_source);
unwarpJImageModel source =
new unwarpJImageModel(sourceImp.getProcessor(), false);
source.setPyramidDepth(0);
source.getThread().start();
// Load transformation
int intervals=unwarpJMiscTools.
numberOfIntervalsOfTransformation(fn_tnf);
double [][]cx=new double[intervals+3][intervals+3];
double [][]cy=new double[intervals+3][intervals+3];
unwarpJMiscTools.loadTransformation(fn_tnf, cx, cy);
// Join threads
try {
source.getThread().join();
} catch (InterruptedException e) {
IJ.error("Unexpected interruption exception " + e);
}
// Apply transformation to source
unwarpJMiscTools.applyTransformationToSource(
sourceImp,targetImp,source,intervals,cx,cy);
// Save results
FileSaver fs=new FileSaver(sourceImp);
if (jpeg_output) {
JpegWriter js = new JpegWriter();
js.setQuality(100);
WindowManager.setTempCurrentImage(sourceImp);
js.run(fn_out);
} else
fs.saveAsTiff(fn_out);
}
} /* end class UnwarpJ_ */
/*====================================================================
| unwarpJClearAll
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJClearAll
extends
Dialog
implements
ActionListener
{ /* begin class unwarpJClearAll */
/*....................................................................
Private variables
....................................................................*/
private ImagePlus sourceImp;
private ImagePlus targetImp;
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void actionPerformed (
final ActionEvent ae
) {
if (ae.getActionCommand().equals("Clear All")) {
sourcePh.removePoints();
targetPh.removePoints();
setVisible(false);
}
else if (ae.getActionCommand().equals("Cancel")) {
setVisible(false);
}
} /* end actionPerformed */
/*------------------------------------------------------------------*/
public Insets getInsets (
) {
return(new Insets(0, 20, 20, 20));
} /* end getInsets */
/*------------------------------------------------------------------*/
unwarpJClearAll (
final Frame parentWindow,
final ImagePlus sourceImp,
final ImagePlus targetImp,
final unwarpJPointHandler sourcePh,
final unwarpJPointHandler targetPh
) {
super(parentWindow, "Removing Points", true);
this.sourceImp = sourceImp;
this.targetImp = targetImp;
this.sourcePh = sourcePh;
this.targetPh = targetPh;
setLayout(new GridLayout(0, 1));
final Button removeButton = new Button("Clear All");
removeButton.addActionListener(this);
final Button cancelButton = new Button("Cancel");
cancelButton.addActionListener(this);
final Label separation1 = new Label("");
final Label separation2 = new Label("");
add(separation1);
add(removeButton);
add(separation2);
add(cancelButton);
pack();
} /* end unwarpJClearAll */
} /* end class unwarpJClearAll */
/*====================================================================
| unwarpJCredits
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJCredits
extends
Dialog
{ /* begin class unwarpJCredits */
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public Insets getInsets (
) {
return(new Insets(0, 20, 20, 20));
} /* end getInsets */
/*------------------------------------------------------------------*/
public unwarpJCredits (
final Frame parentWindow
) {
super(parentWindow, "UnwarpJ", true);
setLayout(new BorderLayout(0, 20));
final Label separation = new Label("");
final Panel buttonPanel = new Panel();
buttonPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Button doneButton = new Button("Done");
doneButton.addActionListener(
new ActionListener (
) {
public void actionPerformed (
final ActionEvent ae
) {
if (ae.getActionCommand().equals("Done")) {
dispose();
}
}
}
);
buttonPanel.add(doneButton);
final TextArea text = new TextArea(22, 72);
text.setEditable(false);
text.append("\n");
text.append(" This work is based on the following paper:\n");
text.append("\n");
text.append(" C.O.S. Sorzano, P. Th" + (char)233 + "venaz, M. Unser\n");
text.append(" Elastic Registration of Biological Images Using Vector-Spline Regularization\n");
text.append(" IEEE Transactions on Biomedical Engineering\n");
text.append(" vol. ??, no. ??, pp. ??-??, July 2005.\n");
text.append("\n");
text.append(" This paper is available on-line at\n");
text.append(" http://bigwww.epfl.ch/publications/sorzano0501.html\n");
text.append("\n");
text.append(" Other relevant on-line publications are available at\n");
text.append(" http://bigwww.epfl.ch/publications/\n");
text.append("\n");
text.append(" Additional help available at\n");
text.append(" http://bigwww.epfl.ch/thevenaz/UnwarpJ/\n");
text.append("\n");
text.append(" You'll be free to use this software for research purposes, but\n");
text.append(" you should not redistribute it without our consent. In addition,\n");
text.append(" we expect you to include a citation or acknowledgment whenever\n");
text.append(" you present or publish results that are based on it.\n");
add("North", separation);
add("Center", text);
add("South", buttonPanel);
pack();
} /* end unwarpJCredits */
} /* end class unwarpJCredits */
/*====================================================================
| unwarpJCumulativeQueue
\===================================================================*/
class unwarpJCumulativeQueue extends Vector {
private int ridx;
private int widx;
private int currentLength;
private double sum;
/*------------------------------------------------------------------*/
public unwarpJCumulativeQueue(int length)
{currentLength=ridx=widx=0; setSize(length);}
/*------------------------------------------------------------------*/
public int currentSize(){return currentLength;}
/*------------------------------------------------------------------*/
public double getSum(){return sum;}
/*------------------------------------------------------------------*/
public double pop_front() {
if (currentLength==0) return 0.0;
double x=((Double)elementAt(ridx)).doubleValue();
currentLength--;
sum-=x;
ridx++; if (ridx==size()) ridx=0;
return x;
}
/*------------------------------------------------------------------*/
public void push_back(double x) {
if (currentLength==size()) pop_front();
setElementAt(new Double(x),widx);
currentLength++;
sum+=x;
widx++; if (widx==size()) widx=0;
}
} /* end class unwarpJCumulativeQueue */
/*====================================================================
| unwarpJDialog
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJDialog
extends
Dialog
implements
ActionListener
{ /* begin class unwarpJDialog */
/*....................................................................
Private variables
....................................................................*/
// Advanced dialog
private Dialog advanced_dlg = null;
// List of available images in ImageJ
private ImagePlus[] imageList;
// Image representations (canvas and ImagePlus)
private ImageCanvas sourceIc;
private ImageCanvas targetIc;
private ImagePlus sourceImp;
private ImagePlus targetImp;
// Image models
private unwarpJImageModel source;
private unwarpJImageModel target;
// Image Masks
private unwarpJMask sourceMsk;
private unwarpJMask targetMsk;
// Point handlers for the landmarks
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
// Toolbar handler
private boolean clearMask=false;
private unwarpJPointToolbar tb
= new unwarpJPointToolbar(Toolbar.getInstance(),this);
// Final action
private boolean finalActionLaunched=false;
private boolean stopRegistration=false;
// Dialog related
private final Button DoneButton = new Button("Done");
private TextField min_scaleDeformationTextField;
private TextField max_scaleDeformationTextField;
private TextField divWeightTextField;
private TextField curlWeightTextField;
private TextField landmarkWeightTextField;
private TextField imageWeightTextField;
private TextField stopThresholdTextField;
private int sourceChoiceIndex = 0;
private int targetChoiceIndex = 1;
private static int min_scale_deformation = 0;
private static int max_scale_deformation = 2;
private static int mode = 1;
private Checkbox ckRichOutput;
private Checkbox ckSaveTransformation;
// Transformation parameters
private static int MIN_SIZE = 8;
private static double divWeight = 0;
private static double curlWeight = 0;
private static double landmarkWeight = 0;
private static double imageWeight = 1;
private static boolean richOutput = false;
private static boolean saveTransformation = false;
private static int min_scale_image = 0;
private static int imagePyramidDepth = 3;
private static double stopThreshold = 1e-2;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void actionPerformed (
final ActionEvent ae
) {
if (ae.getActionCommand().equals("Cancel")) {
dispose();
restoreAll();
}
else if (ae.getActionCommand().equals("Done")) {
dispose();
joinThreads();
imagePyramidDepth = max_scale_deformation-min_scale_deformation+1;
divWeight = Double.valueOf(divWeightTextField.getText()).doubleValue();
curlWeight = Double.valueOf(curlWeightTextField.getText()).doubleValue();
landmarkWeight = Double.valueOf(landmarkWeightTextField.getText()).doubleValue();
imageWeight = Double.valueOf(imageWeightTextField.getText()).doubleValue();
//richOutput = ckRichOutput.getState();
richOutput = false;
saveTransformation = ckSaveTransformation.getState();
int outputLevel=1;
boolean showMarquardtOptim=false;
if (richOutput) {
outputLevel++;
showMarquardtOptim=true;
}
unwarpJFinalAction finalAction =
new unwarpJFinalAction(this);
finalAction.setup(sourceImp,targetImp,
source, target, sourcePh, targetPh,
sourceMsk, targetMsk, min_scale_deformation, max_scale_deformation,
min_scale_image, divWeight, curlWeight, landmarkWeight, imageWeight,
stopThreshold, outputLevel, showMarquardtOptim, mode);
finalActionLaunched=true;
tb.setAllUp(); tb.repaint();
finalAction.getThread().start();
}
else if (ae.getActionCommand().equals("Credits...")) {
final unwarpJCredits dialog = new unwarpJCredits(IJ.getInstance());
GUI.center(dialog);
dialog.setVisible(true);
} else if (ae.getActionCommand().equals("Advanced Options")) {
advanced_dlg.setVisible(true);
} else if (ae.getActionCommand().equals("Done")) {
advanced_dlg.setVisible(false);
}
} /* end actionPerformed */
/*------------------------------------------------------------------*/
public void applyTransformationToSource(
int intervals,
double [][]cx,
double [][]cy) {
// Apply transformation
unwarpJMiscTools.applyTransformationToSource(
sourceImp,targetImp,source,intervals,cx,cy);
// Restart the computation of the model
cancelSource();
targetPh.removePoints();
createSourceImage();
setSecondaryPointHandlers();
}
/*------------------------------------------------------------------*/
public void createAdvancedOptions() {
advanced_dlg = new Dialog(new Frame(), "Advanced Options", true);
// Create min_scale_deformation, max_scale_deformation panel
advanced_dlg.setLayout(new GridLayout(0, 1));
final Choice min_scale_deformationChoice = new Choice();
final Choice max_scale_deformationChoice = new Choice();
final Panel min_scale_deformationPanel = new Panel();
final Panel max_scale_deformationPanel = new Panel();
min_scale_deformationPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
max_scale_deformationPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label min_scale_deformationLabel = new Label("Initial Deformation: ");
final Label max_scale_deformationLabel = new Label("Final Deformation: ");
min_scale_deformationChoice.add("Very Coarse");
min_scale_deformationChoice.add("Coarse");
min_scale_deformationChoice.add("Fine");
min_scale_deformationChoice.add("Very Fine");
max_scale_deformationChoice.add("Very Coarse");
max_scale_deformationChoice.add("Coarse");
max_scale_deformationChoice.add("Fine");
max_scale_deformationChoice.add("Very Fine");
min_scale_deformationChoice.select(min_scale_deformation);
max_scale_deformationChoice.select(max_scale_deformation);
min_scale_deformationChoice.addItemListener(
new ItemListener (
) {
public void itemStateChanged (
final ItemEvent ie
) {
final int new_min_scale_deformation =
min_scale_deformationChoice.getSelectedIndex();
int new_max_scale_deformation=max_scale_deformation;
if (max_scale_deformation<new_min_scale_deformation)
new_max_scale_deformation=new_min_scale_deformation;
if (new_min_scale_deformation!=min_scale_deformation ||
new_max_scale_deformation!=max_scale_deformation) {
min_scale_deformation=new_min_scale_deformation;
max_scale_deformation=new_max_scale_deformation;
computeImagePyramidDepth();
restartModelThreads();
}
min_scale_deformationChoice.select(min_scale_deformation);
max_scale_deformationChoice.select(max_scale_deformation);
}
}
);
max_scale_deformationChoice.addItemListener(
new ItemListener (
) {
public void itemStateChanged (
final ItemEvent ie
) {
final int new_max_scale_deformation =
max_scale_deformationChoice.getSelectedIndex();
int new_min_scale_deformation=min_scale_deformation;
if (new_max_scale_deformation<min_scale_deformation)
new_min_scale_deformation=new_max_scale_deformation;
if (new_max_scale_deformation!=max_scale_deformation ||
new_min_scale_deformation!=min_scale_deformation) {
min_scale_deformation=new_min_scale_deformation;
max_scale_deformation=new_max_scale_deformation;
computeImagePyramidDepth();
restartModelThreads();
}
max_scale_deformationChoice.select(max_scale_deformation);
min_scale_deformationChoice.select(min_scale_deformation);
}
}
);
min_scale_deformationPanel.add(min_scale_deformationLabel);
max_scale_deformationPanel.add(max_scale_deformationLabel);
min_scale_deformationPanel.add(min_scale_deformationChoice);
max_scale_deformationPanel.add(max_scale_deformationChoice);
// Create div and curl weight panels
final Panel divWeightPanel= new Panel();
divWeightPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label label_divWeight = new Label();
label_divWeight.setText("Divergence Weight:");
divWeightTextField = new TextField("", 5);
divWeightTextField.setText(""+divWeight);
divWeightTextField.addTextListener(
new TextListener (
) {
public void textValueChanged (
final TextEvent e
) {
boolean validNumber =true;
try {
divWeight = Double.valueOf(divWeightTextField.getText()).doubleValue();
} catch (NumberFormatException n) {
validNumber = false;
}
DoneButton.setEnabled(validNumber);
}
}
);
divWeightPanel.add(label_divWeight);
divWeightPanel.add(divWeightTextField);
divWeightPanel.setVisible(true);
final Panel curlWeightPanel= new Panel();
curlWeightPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label label_curlWeight = new Label();
label_curlWeight.setText("Curl Weight:");
curlWeightTextField = new TextField("", 5);
curlWeightTextField.setText(""+curlWeight);
curlWeightTextField.addTextListener(
new TextListener (
) {
public void textValueChanged (
final TextEvent e
) {
boolean validNumber =true;
try {
curlWeight = Double.valueOf(curlWeightTextField.getText()).doubleValue();
} catch (NumberFormatException n) {
validNumber = false;
}
DoneButton.setEnabled(validNumber);
}
}
);
curlWeightPanel.add(label_curlWeight);
curlWeightPanel.add(curlWeightTextField);
curlWeightPanel.setVisible(true);
// Create landmark and image weight panels
final Panel landmarkWeightPanel= new Panel();
landmarkWeightPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label label_landmarkWeight = new Label();
label_landmarkWeight.setText("Landmark Weight:");
landmarkWeightTextField = new TextField("", 5);
landmarkWeightTextField.setText(""+landmarkWeight);
landmarkWeightTextField.addTextListener(
new TextListener (
) {
public void textValueChanged (
final TextEvent e
) {
boolean validNumber =true;
try {
landmarkWeight = Double.valueOf(landmarkWeightTextField.getText()).doubleValue();
} catch (NumberFormatException n) {
validNumber = false;
}
DoneButton.setEnabled(validNumber);
}
}
);
landmarkWeightPanel.add(label_landmarkWeight);
landmarkWeightPanel.add(landmarkWeightTextField);
landmarkWeightPanel.setVisible(true);
final Panel imageWeightPanel= new Panel();
imageWeightPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label label_imageWeight = new Label();
label_imageWeight.setText("Image Weight:");
imageWeightTextField = new TextField("", 5);
imageWeightTextField.setText(""+imageWeight);
imageWeightTextField.addTextListener(
new TextListener (
) {
public void textValueChanged (
final TextEvent e
) {
boolean validNumber =true;
try {
imageWeight = Double.valueOf(imageWeightTextField.getText()).doubleValue();
} catch (NumberFormatException n) {
validNumber = false;
}
DoneButton.setEnabled(validNumber);
}
}
);
imageWeightPanel.add(label_imageWeight);
imageWeightPanel.add(imageWeightTextField);
imageWeightPanel.setVisible(true);
// Create stopThreshold panel
final Panel stopThresholdPanel= new Panel();
stopThresholdPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label label_stopThreshold = new Label();
label_stopThreshold.setText("Stop Threshold:");
stopThresholdTextField = new TextField("", 5);
stopThresholdTextField.setText(""+stopThreshold);
stopThresholdTextField.addTextListener(
new TextListener (
) {
public void textValueChanged (
final TextEvent e
) {
boolean validNumber =true;
try {
stopThreshold = Double.valueOf(stopThresholdTextField.getText()).doubleValue();
} catch (NumberFormatException n) {
validNumber = false;
}
DoneButton.setEnabled(validNumber);
}
}
);
stopThresholdPanel.add(label_stopThreshold);
stopThresholdPanel.add(stopThresholdTextField);
stopThresholdPanel.setVisible(true);
// Create checkbox for creating rich output
final Panel richOutputPanel=new Panel();
richOutputPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
ckRichOutput=new Checkbox("Verbose",richOutput);
richOutputPanel.add(ckRichOutput);
// Create checkbox for saving the transformation
final Panel saveTransformationPanel=new Panel();
saveTransformationPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
ckSaveTransformation=new Checkbox("Save Transformation",saveTransformation);
saveTransformationPanel.add(ckSaveTransformation);
// Create close button
final Panel buttonPanel = new Panel();
buttonPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Button CloseButton = new Button("Close");
CloseButton.addActionListener(
new ActionListener (
) {
public void actionPerformed (
final ActionEvent ae
) {
if (ae.getActionCommand().equals("Close")) {
advanced_dlg.dispose();
}
}
}
);
buttonPanel.add(CloseButton);
// Build separations
final Label separation1 = new Label("");
// Create the dialog
advanced_dlg.add(min_scale_deformationPanel);
advanced_dlg.add(max_scale_deformationPanel);
advanced_dlg.add(divWeightPanel);
advanced_dlg.add(curlWeightPanel);
advanced_dlg.add(landmarkWeightPanel);
advanced_dlg.add(imageWeightPanel);
advanced_dlg.add(stopThresholdPanel);
advanced_dlg.add(richOutputPanel);
advanced_dlg.add(saveTransformationPanel);
advanced_dlg.add(separation1);
advanced_dlg.add(buttonPanel);
advanced_dlg.pack();
advanced_dlg.setVisible(false);
}
/*------------------------------------------------------------------*/
public void freeMemory() {
advanced_dlg = null;
imageList = null;
sourceIc = null;
targetIc = null;
sourceImp = null;
targetImp = null;
source = null;
target = null;
sourcePh = null;
targetPh = null;
tb = null;
Runtime.getRuntime().gc();
}
/*------------------------------------------------------------------*/
public void grayImage(final unwarpJPointHandler ph) {
if (ph==sourcePh) {
int Xdim=source.getWidth();
int Ydim=source.getHeight();
FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
int ij=0;
double []source_data=source.getImage();
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++,ij++)
if (sourceMsk.getValue(j,i))
fp.putPixelValue(j,i, source_data[ij]);
else fp.putPixelValue(j,i,0.5*source_data[ij]);
fp.resetMinAndMax();
sourceImp.setProcessor(sourceImp.getTitle(),fp);
sourceImp.updateImage();
} else {
int Xdim=target.getWidth();
int Ydim=target.getHeight();
FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
double []target_data=target.getImage();
int ij=0;
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++,ij++)
if (targetMsk.getValue(j,i))
fp.putPixelValue(j,i, target_data[ij]);
else fp.putPixelValue(j,i,0.5*target_data[ij]);
fp.resetMinAndMax();
targetImp.setProcessor(targetImp.getTitle(),fp);
targetImp.updateImage();
}
}
/*------------------------------------------------------------------*/
public boolean isFinalActionLaunched () {return finalActionLaunched;}
/*------------------------------------------------------------------*/
public boolean isClearMaskSet () {return clearMask;}
/*------------------------------------------------------------------*/
public boolean isSaveTransformationSet () {return saveTransformation;}
/*------------------------------------------------------------------*/
public boolean isStopRegistrationSet () {return stopRegistration;}
/*------------------------------------------------------------------*/
public Insets getInsets (
) {
return(new Insets(0, 20, 20, 20));
} /* end getInsets */
/*------------------------------------------------------------------*/
public void joinThreads (
) {
try {
source.getThread().join();
target.getThread().join();
} catch (InterruptedException e) {
IJ.error("Unexpected interruption exception" + e);
}
} /* end joinSourceThread */
/*------------------------------------------------------------------*/
public void restoreAll (
) {
cancelSource();
cancelTarget();
tb.restorePreviousToolbar();
Toolbar.getInstance().repaint();
unwarpJProgressBar.resetProgressBar();
Runtime.getRuntime().gc();
} /* end restoreAll */
/*------------------------------------------------------------------*/
public void setClearMask (boolean val) {clearMask=val;}
/*------------------------------------------------------------------*/
public void setStopRegistration () {stopRegistration=true;}
/*------------------------------------------------------------------*/
public unwarpJDialog (
final Frame parentWindow,
final ImagePlus[] imageList
) {
super(parentWindow, "UnwarpJ", false);
this.imageList = imageList;
// Start concurrent image processing threads
createSourceImage();
createTargetImage();
setSecondaryPointHandlers();
// Create Source panel
setLayout(new GridLayout(0, 1));
final Choice sourceChoice = new Choice();
final Choice targetChoice = new Choice();
final Panel sourcePanel = new Panel();
sourcePanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label sourceLabel = new Label("Source: ");
addImageList(sourceChoice);
sourceChoice.select(sourceChoiceIndex);
sourceChoice.addItemListener(
new ItemListener (
) {
public void itemStateChanged (
final ItemEvent ie
) {
final int newChoiceIndex = sourceChoice.getSelectedIndex();
if (sourceChoiceIndex != newChoiceIndex) {
stopSourceThread();
if (targetChoiceIndex != newChoiceIndex) {
sourceChoiceIndex = newChoiceIndex;
cancelSource();
targetPh.removePoints();
createSourceImage();
setSecondaryPointHandlers();
}
else {
stopTargetThread();
targetChoiceIndex = sourceChoiceIndex;
sourceChoiceIndex = newChoiceIndex;
targetChoice.select(targetChoiceIndex);
permuteImages();
}
}
repaint();
}
}
);
sourcePanel.add(sourceLabel);
sourcePanel.add(sourceChoice);
// Create target panel
final Panel targetPanel = new Panel();
targetPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label targetLabel = new Label("Target: ");
addImageList(targetChoice);
targetChoice.select(targetChoiceIndex);
targetChoice.addItemListener(
new ItemListener (
) {
public void itemStateChanged (
final ItemEvent ie
) {
final int newChoiceIndex = targetChoice.getSelectedIndex();
if (targetChoiceIndex != newChoiceIndex) {
stopTargetThread();
if (sourceChoiceIndex != newChoiceIndex) {
targetChoiceIndex = newChoiceIndex;
cancelTarget();
sourcePh.removePoints();
createTargetImage();
setSecondaryPointHandlers();
}
else {
stopSourceThread();
sourceChoiceIndex = targetChoiceIndex;
targetChoiceIndex = newChoiceIndex;
sourceChoice.select(sourceChoiceIndex);
permuteImages();
}
}
repaint();
}
}
);
targetPanel.add(targetLabel);
targetPanel.add(targetChoice);
// Create mode panel
setLayout(new GridLayout(0, 1));
final Choice modeChoice = new Choice();
final Panel modePanel = new Panel();
modePanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Label modeLabel = new Label("Registration Mode: ");
modeChoice.add("Fast");
modeChoice.add("Accurate");
modeChoice.select(mode);
modeChoice.addItemListener(
new ItemListener (
) {
public void itemStateChanged (
final ItemEvent ie
) {
final int mode = modeChoice.getSelectedIndex();
if (mode==0) {
// Fast
min_scale_image=1;
} else if (mode==1) {
// Accurate
min_scale_image=0;
}
repaint();
}
}
);
modePanel.add(modeLabel);
modePanel.add(modeChoice);
// Build Advanced Options panel
final Panel advancedOptionsPanel = new Panel();
advancedOptionsPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
final Button advancedOptionsButton = new Button("Advanced Options");
advancedOptionsButton.addActionListener(this);
advancedOptionsPanel.add(advancedOptionsButton);
// Build Done Cancel Credits panel
final Panel buttonPanel = new Panel();
buttonPanel.setLayout(new FlowLayout(FlowLayout.CENTER));
DoneButton.addActionListener(this);
final Button cancelButton = new Button("Cancel");
cancelButton.addActionListener(this);
final Button creditsButton = new Button("Credits...");
creditsButton.addActionListener(this);
buttonPanel.add(cancelButton);
buttonPanel.add(DoneButton);
buttonPanel.add(creditsButton);
// Build separations
final Label separation1 = new Label("");
final Label separation2 = new Label("");
// Finally build dialog
add(separation1);
add(sourcePanel);
add(targetPanel);
add(modePanel);
add(advancedOptionsPanel);
add(separation2);
add(buttonPanel);
pack();
createAdvancedOptions();
} /* end unwarpJDialog */
/*------------------------------------------------------------------*/
public void ungrayImage(final unwarpJPointAction pa) {
if (pa==sourcePh.getPointAction()) {
int Xdim=source.getWidth();
int Ydim=source.getHeight();
FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
int ij=0;
double []source_data=source.getImage();
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++,ij++)
fp.putPixelValue(j,i,source_data[ij]);
fp.resetMinAndMax();
sourceImp.setProcessor(sourceImp.getTitle(),fp);
sourceImp.updateImage();
} else {
int Xdim=target.getWidth();
int Ydim=target.getHeight();
FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
double []target_data=target.getImage();
int ij=0;
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++,ij++)
fp.putPixelValue(j,i,target_data[ij]);
fp.resetMinAndMax();
targetImp.setProcessor(targetImp.getTitle(),fp);
targetImp.updateImage();
}
}
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void addImageList (
final Choice choice
) {
for (int k = 0; (k < imageList.length); k++)
choice.add(imageList[k].getTitle());
} /* end addImageList */
/*------------------------------------------------------------------*/
private void cancelSource (
) {
sourcePh.killListeners();
sourcePh = null;
sourceIc = null;
sourceImp.killRoi();
sourceImp = null;
source = null;
sourceMsk = null;
Runtime.getRuntime().gc();
} /* end cancelSource */
/*------------------------------------------------------------------*/
private void cancelTarget (
) {
targetPh.killListeners();
targetPh = null;
targetIc = null;
targetImp.killRoi();
targetImp = null;
target = null;
targetMsk = null;
Runtime.getRuntime().gc();
} /* end cancelTarget */
/*------------------------------------------------------------------*/
private void computeImagePyramidDepth (
) {
imagePyramidDepth=max_scale_deformation-min_scale_deformation+1;
}
/*------------------------------------------------------------------*/
private void createSourceImage (
) {
sourceImp = imageList[sourceChoiceIndex];
sourceImp.setSlice(1);
source =
new unwarpJImageModel(sourceImp.getProcessor(), false);
source.setPyramidDepth(imagePyramidDepth+min_scale_image);
source.getThread().start();
sourceIc = sourceImp.getWindow().getCanvas();
if (sourceImp.getStackSize()==1) {
// Create an empy mask
sourceMsk = new unwarpJMask(sourceImp.getProcessor(),false);
} else {
// Take the mask from the second slice
sourceImp.setSlice(2);
sourceMsk = new unwarpJMask(sourceImp.getProcessor(),true);
sourceImp.setSlice(1);
}
sourcePh = new unwarpJPointHandler(sourceImp, tb, sourceMsk, this);
tb.setSource(sourceImp, sourcePh);
} /* end createSourceImage */
/*------------------------------------------------------------------*/
private void createTargetImage (
) {
targetImp = imageList[targetChoiceIndex];
targetImp.setSlice(1);
target =
new unwarpJImageModel(targetImp.getProcessor(), true);
target.setPyramidDepth(imagePyramidDepth+min_scale_image);
target.getThread().start();
targetIc = targetImp.getWindow().getCanvas();
if (targetImp.getStackSize()==1) {
// Create an empy mask
targetMsk = new unwarpJMask(targetImp.getProcessor(),false);
} else {
// Take the mask from the second slice
targetImp.setSlice(2);
targetMsk = new unwarpJMask(targetImp.getProcessor(),true);
targetImp.setSlice(1);
}
targetPh = new unwarpJPointHandler(targetImp, tb, targetMsk, this);
tb.setTarget(targetImp, targetPh);
} /* end createTargetImage */
/*------------------------------------------------------------------*/
private void permuteImages (
) {
// Swap image canvas
final ImageCanvas swapIc = sourceIc;
sourceIc = targetIc;
targetIc = swapIc;
// Swap ImagePlus
final ImagePlus swapImp = sourceImp;
sourceImp = targetImp;
targetImp = swapImp;
// Swap Mask
final unwarpJMask swapMsk = sourceMsk;
sourceMsk=targetMsk;
targetMsk=swapMsk;
// Swap Point Handlers
final unwarpJPointHandler swapPh = sourcePh;
sourcePh = targetPh;
targetPh = swapPh;
setSecondaryPointHandlers();
// Inform the Toolbar about the change
tb.setSource(sourceImp, sourcePh);
tb.setTarget(targetImp, targetPh);
// Restart the computation with each image
source =
new unwarpJImageModel(sourceImp.getProcessor(), false);
source.setPyramidDepth(imagePyramidDepth+min_scale_image);
source.getThread().start();
target =
new unwarpJImageModel(targetImp.getProcessor(), true);
target.setPyramidDepth(imagePyramidDepth+min_scale_image);
target.getThread().start();
} /* end permuteImages */
/*------------------------------------------------------------------*/
private void removePoints (
) {
sourcePh.removePoints();
targetPh.removePoints();
}
/*------------------------------------------------------------------*/
private void restartModelThreads (
) {
// Stop threads
stopSourceThread();
stopTargetThread();
// Remove the current image models
source = null;
target = null;
Runtime.getRuntime().gc();
// Now restart the threads
source =
new unwarpJImageModel(sourceImp.getProcessor(), false);
source.setPyramidDepth(imagePyramidDepth+min_scale_image);
source.getThread().start();
target =
new unwarpJImageModel(targetImp.getProcessor(), true);
target.setPyramidDepth(imagePyramidDepth+min_scale_image);
target.getThread().start();
}
/*------------------------------------------------------------------*/
private void setSecondaryPointHandlers (
) {
sourcePh.setSecondaryPointHandler(targetImp, targetPh);
targetPh.setSecondaryPointHandler(sourceImp, sourcePh);
} /* end setSecondaryPointHandler */
/*------------------------------------------------------------------*/
private void stopSourceThread (
) {
// Stop the source image model
while (source.getThread().isAlive()) {
source.getThread().interrupt();
}
source.getThread().interrupted();
} /* end stopSourceThread */
/*------------------------------------------------------------------*/
private void stopTargetThread (
) {
// Stop the target image model
while (target.getThread().isAlive()) {
target.getThread().interrupt();
}
target.getThread().interrupted();
} /* end stopTargetThread */
} /* end class unwarpJDialog */
/*====================================================================
| unwarpJFile
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJFile
extends
Dialog
implements
ActionListener
{ /* begin class unwarpJFile */
/*....................................................................
Private variables
....................................................................*/
private final CheckboxGroup choice = new CheckboxGroup();
private ImagePlus sourceImp;
private ImagePlus targetImp;
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
private unwarpJDialog dialog;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void actionPerformed (
final ActionEvent ae
) {
this.setVisible(false);
if (ae.getActionCommand().equals("Save Landmarks As...")) {
savePoints();
}
else if (ae.getActionCommand().equals("Load Landmarks...")) {
loadPoints();
}
else if (ae.getActionCommand().equals("Show Landmarks")) {
showPoints();
}
else if (ae.getActionCommand().equals("Load Transformation")) {
loadTransformation();
}
else if (ae.getActionCommand().equals("Cancel")) {
}
} /* end actionPerformed */
/*------------------------------------------------------------------*/
public Insets getInsets (
) {
return(new Insets(0, 20, 20, 20));
} /* end getInsets */
/*------------------------------------------------------------------*/
unwarpJFile (
final Frame parentWindow,
final ImagePlus sourceImp,
final ImagePlus targetImp,
final unwarpJPointHandler sourcePh,
final unwarpJPointHandler targetPh,
final unwarpJDialog dialog
) {
super(parentWindow, "I/O Menu", true);
this.sourceImp = sourceImp;
this.targetImp = targetImp;
this.sourcePh = sourcePh;
this.targetPh = targetPh;
this.dialog = dialog;
setLayout(new GridLayout(0, 1));
final Button saveAsButton = new Button("Save Landmarks As...");
final Button loadButton = new Button("Load Landmarks...");
final Button show_PointsButton = new Button("Show Landmarks");
final Button loadTransfButton = new Button("Load Transformation");
final Button cancelButton = new Button("Cancel");
saveAsButton.addActionListener(this);
loadButton.addActionListener(this);
show_PointsButton.addActionListener(this);
loadTransfButton.addActionListener(this);
cancelButton.addActionListener(this);
final Label separation1 = new Label("");
final Label separation2 = new Label("");
add(separation1);
add(loadButton);
add(saveAsButton);
add(show_PointsButton);
add(loadTransfButton);
add(separation2);
add(cancelButton);
pack();
} /* end unwarpJFile */
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void loadPoints (
) {
final Frame f = new Frame();
final FileDialog fd = new FileDialog(f, "Load Points", FileDialog.LOAD);
fd.setVisible(true);
final String path = fd.getDirectory();
final String filename = fd.getFile();
if ((path == null) || (filename == null)) return;
Stack sourceStack = new Stack();
Stack targetStack = new Stack();
unwarpJMiscTools.loadPoints(path+filename,sourceStack,targetStack);
sourcePh.removePoints();
targetPh.removePoints();
while ((!sourceStack.empty()) && (!targetStack.empty())) {
Point sourcePoint = (Point)sourceStack.pop();
Point targetPoint = (Point)targetStack.pop();
sourcePh.addPoint(sourcePoint.x, sourcePoint.y);
targetPh.addPoint(targetPoint.x, targetPoint.y);
}
} /* end loadPoints */
/*------------------------------------------------------------------*/
private void loadTransformation (
) {
final Frame f = new Frame();
final FileDialog fd = new FileDialog(f, "Load Transformation", FileDialog.LOAD);
fd.setVisible(true);
final String path = fd.getDirectory();
final String filename = fd.getFile();
if ((path == null) || (filename == null)) {
return;
}
String fn_tnf=path+filename;
int intervals=unwarpJMiscTools.
numberOfIntervalsOfTransformation(fn_tnf);
double [][]cx=new double[intervals+3][intervals+3];
double [][]cy=new double[intervals+3][intervals+3];
unwarpJMiscTools.loadTransformation(fn_tnf, cx, cy);
// Apply transformation
dialog.applyTransformationToSource(intervals,cx,cy);
}
/*------------------------------------------------------------------*/
private void savePoints (
) {
final Frame f = new Frame();
final FileDialog fd = new FileDialog(f, "Save Points", FileDialog.SAVE);
String filename = targetImp.getTitle();
int dot = filename.lastIndexOf('.');
if (dot == -1) {
fd.setFile(filename + ".txt");
}
else {
filename = filename.substring(0, dot);
fd.setFile(filename + ".txt");
}
fd.setVisible(true);
final String path = fd.getDirectory();
filename = fd.getFile();
if ((path == null) || (filename == null)) {
return;
}
try {
final FileWriter fw = new FileWriter(path + filename);
final Vector sourceList = sourcePh.getPoints();
final Vector targetList = targetPh.getPoints();
Point sourcePoint;
Point targetPoint;
String n;
String xSource;
String ySource;
String xTarget;
String yTarget;
fw.write("Index\txSource\tySource\txTarget\tyTarget\n");
for (int k = 0; (k < sourceList.size()); k++) {
n = "" + k;
while (n.length() < 5) {
n = " " + n;
}
sourcePoint = (Point)sourceList.elementAt(k);
xSource = "" + sourcePoint.x;
while (xSource.length() < 7) {
xSource = " " + xSource;
}
ySource = "" + sourcePoint.y;
while (ySource.length() < 7) {
ySource = " " + ySource;
}
targetPoint = (Point)targetList.elementAt(k);
xTarget = "" + targetPoint.x;
while (xTarget.length() < 7) {
xTarget = " " + xTarget;
}
yTarget = "" + targetPoint.y;
while (yTarget.length() < 7) {
yTarget = " " + yTarget;
}
fw.write(n + "\t" + xSource + "\t" + ySource + "\t" + xTarget + "\t" + yTarget + "\n");
}
fw.close();
} catch (IOException e) {
IJ.error("IOException exception" + e);
} catch (SecurityException e) {
IJ.error("Security exception" + e);
}
} /* end savePoints */
/*------------------------------------------------------------------*/
private void showPoints (
) {
final Vector sourceList = sourcePh.getPoints();
final Vector targetList = targetPh.getPoints();
Point sourcePoint;
Point targetPoint;
String n;
String xTarget;
String yTarget;
String xSource;
String ySource;
IJ.getTextPanel().setFont(new Font("Monospaced", Font.PLAIN, 12));
IJ.setColumnHeadings("Index\txSource\tySource\txTarget\tyTarget");
for (int k = 0; (k < sourceList.size()); k++) {
n = "" + k;
while (n.length() < 5) {
n = " " + n;
}
sourcePoint = (Point)sourceList.elementAt(k);
xTarget = "" + sourcePoint.x;
while (xTarget.length() < 7) {
xTarget = " " + xTarget;
}
yTarget = "" + sourcePoint.y;
while (yTarget.length() < 7) {
yTarget = " " + yTarget;
}
targetPoint = (Point)targetList.elementAt(k);
xSource = "" + targetPoint.x;
while (xSource.length() < 7) {
xSource = " " + xSource;
}
ySource = "" + targetPoint.y;
while (ySource.length() < 7) {
ySource = " " + ySource;
}
IJ.write(n + "\t" + xSource + "\t" + ySource + "\t" + xTarget + "\t" + yTarget);
}
} /* end showPoints */
} /* end class unwarpJFile */
/*====================================================================
| unwarpJFinalAction
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJFinalAction
implements
Runnable
{
/*....................................................................
Private variables
....................................................................*/
private Thread t;
private unwarpJDialog dialog;
// Images
private ImagePlus sourceImp;
private ImagePlus targetImp;
private unwarpJImageModel source;
private unwarpJImageModel target;
// Landmarks
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
// Masks for the images
private unwarpJMask sourceMsk;
private unwarpJMask targetMsk;
// Transformation parameters
private int min_scale_deformation;
private int max_scale_deformation;
private int min_scale_image;
private int outputLevel;
private boolean showMarquardtOptim;
private double divWeight;
private double curlWeight;
private double landmarkWeight;
private double imageWeight;
private double stopThreshold;
private int accurate_mode;
/*....................................................................
Public methods
....................................................................*/
/*********************************************************************
* Return the thread associated with this <code>unwarpJFinalAction</code>
* object.
********************************************************************/
public Thread getThread (
) {
return(t);
} /* end getThread */
/*********************************************************************
* Perform the registration
********************************************************************/
public void run (
) {
// Create output image
int Ydimt=target.getHeight();
int Xdimt=target.getWidth();
int Xdims=source.getWidth();
final FloatProcessor fp=new FloatProcessor(Xdimt,Ydimt);
for (int i=0; i<Ydimt; i++)
for (int j=0; j<Xdimt; j++)
if (sourceMsk.getValue(j,i) && targetMsk.getValue(j,i))
fp.putPixelValue(j,i,(target.getImage())[i*Xdimt+j]-
(source.getImage())[i*Xdims+j]);
else fp.putPixelValue(j,i,0);
fp.resetMinAndMax();
final ImagePlus ip=new ImagePlus("Output", fp);
final ImageWindow iw=new ImageWindow(ip);
ip.updateImage();
// Perform the registration
final unwarpJTransformation warp = new unwarpJTransformation(
sourceImp, targetImp, source, target, sourcePh, targetPh,
sourceMsk, targetMsk, min_scale_deformation, max_scale_deformation,
min_scale_image, divWeight, curlWeight, landmarkWeight, imageWeight,
stopThreshold, outputLevel, showMarquardtOptim, accurate_mode,
dialog.isSaveTransformationSet(), "", ip, dialog);
warp.doRegistration();
dialog.ungrayImage(sourcePh.getPointAction());
dialog.ungrayImage(targetPh.getPointAction());
dialog.restoreAll();
dialog.freeMemory();
}
/*********************************************************************
* Pass parameter from <code>unwarpJDialog</code> to
* <code>unwarpJFinalAction</code>.
********************************************************************/
public void setup (
final ImagePlus sourceImp,
final ImagePlus targetImp,
final unwarpJImageModel source,
final unwarpJImageModel target,
final unwarpJPointHandler sourcePh,
final unwarpJPointHandler targetPh,
final unwarpJMask sourceMsk,
final unwarpJMask targetMsk,
final int min_scale_deformation,
final int max_scale_deformation,
final int min_scale_image,
final double divWeight,
final double curlWeight,
final double landmarkWeight,
final double imageWeight,
final double stopThreshold,
final int outputLevel,
final boolean showMarquardtOptim,
final int accurate_mode
) {
this.sourceImp = sourceImp;
this.targetImp = targetImp;
this.source = source;
this.target = target;
this.sourcePh = sourcePh;
this.targetPh = targetPh;
this.sourceMsk = sourceMsk;
this.targetMsk = targetMsk;
this.min_scale_deformation = min_scale_deformation;
this.max_scale_deformation = max_scale_deformation;
this.min_scale_image = min_scale_image;
this.divWeight = divWeight;
this.curlWeight = curlWeight;
this.landmarkWeight = landmarkWeight;
this.imageWeight = imageWeight;
this.stopThreshold = stopThreshold;
this.outputLevel = outputLevel;
this.showMarquardtOptim = showMarquardtOptim;
this.accurate_mode = accurate_mode;
} /* end setup */
/*********************************************************************
* Start a thread under the control of the main event loop. This thread
* has access to the progress bar, while methods called directly from
* within <code>unwarpJDialog</code> do not because they are
* under the control of its own event loop.
********************************************************************/
public unwarpJFinalAction (
final unwarpJDialog dialog
) {
this.dialog = dialog;
t = new Thread(this);
t.setDaemon(true);
}
}
/*====================================================================
| unwarpJImageModel
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJImageModel
implements
Runnable
{ /* begin class unwarpJImageModel */
// Some constants
private static int min_image_size=4;
/*....................................................................
Private variables
....................................................................*/
// Thread
private Thread t;
// Stack for the pyramid of images/coefficients
private final Stack cpyramid = new Stack();
private final Stack imgpyramid = new Stack();
// Original image, image spline coefficients, and gradient
private double[] image;
private double[] coefficient;
// Current image (the size might be different from the original)
private double[] currentImage;
private double[] currentCoefficient;
private int currentWidth;
private int currentHeight;
private int twiceCurrentWidth;
private int twiceCurrentHeight;
// Size and other information
private int width;
private int height;
private int twiceWidth;
private int twiceHeight;
private int pyramidDepth;
private int currentDepth;
private int smallestWidth;
private int smallestHeight;
private boolean isTarget;
private boolean coefficientsAreMirrored;
// Some variables to speedup interpolation
// All these information is set through prepareForInterpolation()
private double x; // Point to interpolate
private double y;
public int xIndex[]; // Indexes related
public int yIndex[];
private double xWeight[]; // Weights of the splines related
private double yWeight[];
private double dxWeight[]; // Weights of the derivatives splines related
private double dyWeight[];
private double d2xWeight[]; // Weights of the second derivatives splines related
private double d2yWeight[];
private boolean fromCurrent; // Interpolation source (current or original)
private int widthToUse; // Size of the image used for the interpolation
private int heightToUse;
// Some variables to speedup interpolation (precomputed)
// All these information is set through prepareForInterpolation()
public int prec_xIndex[][]; // Indexes related
public int prec_yIndex[][];
private double prec_xWeight[][]; // Weights of the splines related
private double prec_yWeight[][];
private double prec_dxWeight[][]; // Weights of the derivatives splines related
private double prec_dyWeight[][];
private double prec_d2xWeight[][]; // Weights of the second derivatives splines related
private double prec_d2yWeight[][];
/*....................................................................
Public methods
....................................................................*/
/* Clear the pyramid. */
public void clearPyramids (
) {
cpyramid.removeAllElements();
imgpyramid.removeAllElements();
} /* end clearPyramid */
/*------------------------------------------------------------------*/
/* Return the full-size B-spline coefficients. */
public double[] getCoefficient () {return coefficient;}
/*------------------------------------------------------------------*/
/* Return the current height of the image/coefficients. */
public int getCurrentHeight() {return currentHeight;}
/*------------------------------------------------------------------*/
/* Return the current image of the image/coefficients. */
public double[] getCurrentImage() {return currentImage;}
/*------------------------------------------------------------------*/
/* Return the current width of the image/coefficients. */
public int getCurrentWidth () {return currentWidth;}
/*------------------------------------------------------------------*/
/* Return the relationship between the current size of the image
and the original size. */
public double getFactorHeight () {return (double)currentHeight/height;}
/*------------------------------------------------------------------*/
/* Return the relationship between the current size of the image
and the original size. */
public double getFactorWidth () {return (double)currentWidth/width;}
/*------------------------------------------------------------------*/
/* Return the current depth of the image/coefficients. */
public int getCurrentDepth() {return currentDepth;}
/*------------------------------------------------------------------*/
/* Return the full-size image height. */
public int getHeight () {return(height);}
/*------------------------------------------------------------------*/
/* Return the full-size image. */
public double[] getImage () {return image;}
/*------------------------------------------------------------------*/
public double getPixelValFromPyramid(
int x, // Pixel location
int y) {
return currentImage[y*currentWidth+x];
}
/*------------------------------------------------------------------*/
/* Return the depth of the image pyramid. A depth 1 means
that one coarse resolution level is present in the stack. The
full-size level is not placed on the stack. */
public int getPyramidDepth () {return(pyramidDepth);}
/*------------------------------------------------------------------*/
/* Return the height of the smallest image in the pyramid. */
public int getSmallestHeight () {return(smallestHeight);}
/*------------------------------------------------------------------*/
/* Return the width of the smallest image in the pyramid. */
public int getSmallestWidth () {return(smallestWidth);}
/*------------------------------------------------------------------*/
/* Return the thread associated. */
public Thread getThread () {return(t);}
/*------------------------------------------------------------------*/
/* Return the full-size image width. */
public int getWidth () {return(width);}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightDx(int l, int m) {return yWeight[l]*dxWeight[m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightDxDx(int l, int m) {return yWeight[l]*d2xWeight[m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightDxDy(int l, int m) {return dyWeight[l]*dxWeight[m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightDy(int l, int m) {return dyWeight[l]*xWeight[m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightDyDy(int l, int m) {return d2yWeight[l]*xWeight[m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double getWeightI(int l, int m) {return yWeight[l]*xWeight[m];}
/*------------------------------------------------------------------*/
/*
There are two types of interpolation routines. Those that use
precomputed weights and those that don't.
An example of use of the ones without precomputation is the
following:
// Set of B-spline coefficients
double [][]c;
// Set these coefficients to an interpolator
unwarpJImageModel sw = new unwarpJImageModel(c);
// Compute the transformation mapping
for (int v=0; v<ImageHeight; v++) {
final double tv = (double)(v * intervals) / (double)(ImageHeight - 1) + 1.0F;
for (int u = 0; u<ImageeWidth; u++) {
final double tu = (double)(u * intervals) / (double)(ImageWidth - 1) + 1.0F;
sw.prepareForInterpolation(tu, tv, ORIGINAL);
interpolated_val[v][u] = sw.interpolateI();
*/
/*------------------------------------------------------------------*/
/*------------------------------------------------------------------*/
/* Interpolate the X and Y derivatives of the image at a
given point. */
public void interpolateD(double []D) {
// Only SplineDegree=3 is implemented
D[0]=D[1]=0.0F;
for (int j = 0; j<4; j++) {
double sx=0.0F, sy=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1) {
double c;
if (fromCurrent) c=currentCoefficient[p + ix];
else c=coefficient[p + ix];
sx += dxWeight[i]*c;
sy += xWeight[i]*c;
}
}
D[0]+= yWeight[j] * sx;
D[1]+=dyWeight[j] * sy;
}
}
} /* end Interpolate D */
/*------------------------------------------------------------------*/
/* Interpolate the XY, XX and YY derivatives of the image at a
given point. */
public void interpolateD2 (double []D2) {
// Only SplineDegree=3 is implemented
D2[0]=D2[1]=D2[2]=0.0F;
for (int j = 0; j<4; j++) {
double sxy=0.0F, sxx=0.0F, syy=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1) {
double c;
if (fromCurrent) c=currentCoefficient[p + ix];
else c=coefficient[p + ix];
sxy += dxWeight[i]*c;
sxx += d2xWeight[i]*c;
syy += xWeight[i]*c;
}
}
D2[0]+= dyWeight[j] * sxy;
D2[1]+= yWeight[j] * sxx;
D2[2]+=d2yWeight[j] * syy;
}
}
} /* end Interpolate dxdy, dxdx and dydy */
/*------------------------------------------------------------------*/
/* Interpolate the X derivative of the image at a given point. */
public double interpolateDx () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += dxWeight[i]*currentCoefficient[p + ix];
else s += dxWeight[i]*coefficient[p + ix];
}
ival+=yWeight[j] * s;
}
}
return ival;
} /* end Interpolate Dx */
/*------------------------------------------------------------------*/
/* Interpolate the X derivative of the image at a given point. */
public double interpolateDxDx () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += d2xWeight[i]*currentCoefficient[p + ix];
else s += d2xWeight[i]*coefficient[p + ix];
}
ival+=yWeight[j] * s;
}
}
return ival;
} /* end Interpolate DxDx */
/*------------------------------------------------------------------*/
/* Interpolate the X derivative of the image at a given point. */
public double interpolateDxDy () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += dxWeight[i]*currentCoefficient[p + ix];
else s += dxWeight[i]*coefficient[p + ix];
}
ival+=dyWeight[j] * s;
}
}
return ival;
} /* end Interpolate DxDy */
/*------------------------------------------------------------------*/
/* Interpolate the Y derivative of the image at a given point. */
public double interpolateDy () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += xWeight[i]*currentCoefficient[p + ix];
else s += xWeight[i]*coefficient[p + ix];
}
ival+=dyWeight[j] * s;
}
}
return ival;
} /* end Interpolate Dy */
/*------------------------------------------------------------------*/
/* Interpolate the X derivative of the image at a given point. */
public double interpolateDyDy () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += xWeight[i]*currentCoefficient[p + ix];
else s += xWeight[i]*coefficient[p + ix];
}
ival+=d2yWeight[j] * s;
}
}
return ival;
} /* end Interpolate DyDy */
/*------------------------------------------------------------------*/
/* Interpolate the image at a given point. */
public double interpolateI () {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=yIndex[j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=xIndex[i];
if (ix!=-1)
if (fromCurrent) s += xWeight[i]*currentCoefficient[p + ix];
else s += xWeight[i]*coefficient[p + ix];
}
ival+=yWeight[j] * s;
}
}
return ival;
} /* end Interpolate Image */
/*------------------------------------------------------------------*/
public boolean isFinest() {return cpyramid.isEmpty();}
/*------------------------------------------------------------------*/
public void popFromPyramid(
){
// Pop coefficients
if (cpyramid.isEmpty()) {
currentWidth = width;
currentHeight = height;
currentCoefficient = coefficient;
} else {
currentWidth = ((Integer)cpyramid.pop()).intValue();
currentHeight = ((Integer)cpyramid.pop()).intValue();
currentCoefficient = (double [])cpyramid.pop();
}
twiceCurrentWidth = 2*currentWidth;
twiceCurrentHeight = 2*currentHeight;
if (currentDepth>0) currentDepth--;
// Pop image
if (isTarget && !imgpyramid.isEmpty()) {
if (currentWidth != ((Integer)imgpyramid.pop()).intValue())
System.out.println("I cannot understand");
if (currentHeight != ((Integer)imgpyramid.pop()).intValue())
System.out.println("I cannot understand");
currentImage = (double [])imgpyramid.pop();
} else currentImage = image;
}
/*------------------------------------------------------------------*/
/* fromCurrent=true --> The interpolation is prepared to be done
from the current image in the pyramid.
fromCurrent=false --> The interpolation is prepared to be done
from the original image. */
public void prepareForInterpolation(
double x,
double y,
boolean fromCurrent
) {
// Remind this point for interpolation
this.x=x;
this.y=y;
this.fromCurrent=fromCurrent;
if (fromCurrent) {widthToUse=currentWidth; heightToUse=currentHeight;}
else {widthToUse=width; heightToUse=height;}
int ix=(int)x;
int iy=(int)y;
int twiceWidthToUse =2*widthToUse;
int twiceHeightToUse=2*heightToUse;
// Set X indexes
// p is the index of the rightmost influencing spline
int p = (0.0 <= x) ? (ix + 2) : (ix + 1);
for (int k = 0; k<4; p--, k++) {
if (coefficientsAreMirrored) {
int q = (p < 0) ? (-1 - p) : (p);
if (twiceWidthToUse <= q) q -= twiceWidthToUse * (q / twiceWidthToUse);
xIndex[k] = (widthToUse <= q) ? (twiceWidthToUse - 1 - q) : (q);
} else
xIndex[k] = (p<0 || p>=widthToUse) ? (-1):(p);
}
// Set Y indexes
p = (0.0 <= y) ? (iy + 2) : (iy + 1);
for (int k = 0; k<4; p--, k++) {
if (coefficientsAreMirrored) {
int q = (p < 0) ? (-1 - p) : (p);
if (twiceHeightToUse <= q) q -= twiceHeightToUse * (q / twiceHeightToUse);
yIndex[k] = (heightToUse <= q) ? (twiceHeightToUse - 1 - q) : (q);
} else
yIndex[k] = (p<0 || p>=heightToUse) ? (-1):(p);
}
// Compute how much the sample depart from an integer position
double ex = x - ((0.0 <= x) ? (ix) : (ix - 1));
double ey = y - ((0.0 <= y) ? (iy) : (iy - 1));
// Set X weights for the image and derivative interpolation
double s = 1.0F - ex;
dxWeight[0] = 0.5F * ex * ex;
xWeight[0] = ex * dxWeight[0] / 3.0F; // Bspline03(x-ix-2)
dxWeight[3] = -0.5F * s * s;
xWeight[3] = s * dxWeight[3] / -3.0F; // Bspline03(x-ix+1)
dxWeight[1] = 1.0F - 2.0F * dxWeight[0] + dxWeight[3];
//xWeight[1] = 2.0F / 3.0F + (1.0F + ex) * dxWeight[3]; // Bspline03(x-ix-1);
xWeight[1] = unwarpJMathTools.Bspline03(x-ix-1);
dxWeight[2] = 1.5F * ex * (ex - 4.0F/ 3.0F);
xWeight[2] = 2.0F / 3.0F - (2.0F - ex) * dxWeight[0]; // Bspline03(x-ix)
d2xWeight[0] = ex;
d2xWeight[1] = s-2*ex;
d2xWeight[2] = ex-2*s;
d2xWeight[3] = s;
// Set Y weights for the image and derivative interpolation
double t = 1.0F - ey;
dyWeight[0] = 0.5F * ey * ey;
yWeight[0] = ey * dyWeight[0] / 3.0F;
dyWeight[3] = -0.5F * t * t;
yWeight[3] = t * dyWeight[3] / -3.0F;
dyWeight[1] = 1.0F - 2.0F * dyWeight[0] + dyWeight[3];
yWeight[1] = 2.0F / 3.0F + (1.0F + ey) * dyWeight[3];
dyWeight[2] = 1.5F * ey * (ey - 4.0F/ 3.0F);
yWeight[2] = 2.0F / 3.0F - (2.0F - ey) * dyWeight[0];
d2yWeight[0] = ey;
d2yWeight[1] = t-2*ey;
d2yWeight[2] = ey-2*t;
d2yWeight[3] = t;
} /* prepareForInterpolation */
/*------------------------------------------------------------------*/
/* Return the width of the precomputed vectors */
public int precomputed_getWidth() {return prec_yWeight.length;}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightDx(int l, int m, int u, int v)
{return prec_yWeight[v][l]*prec_dxWeight[u][m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightDxDx(int l, int m, int u, int v)
{return prec_yWeight[v][l]*prec_d2xWeight[u][m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightDxDy(int l, int m, int u, int v)
{return prec_dyWeight[v][l]*prec_dxWeight[u][m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightDy(int l, int m, int u, int v)
{return prec_dyWeight[v][l]*prec_xWeight[u][m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightDyDy(int l, int m, int u, int v)
{return prec_d2yWeight[v][l]*prec_xWeight[u][m];}
/*------------------------------------------------------------------*/
/* Return the weight of the coefficient l,m (yIndex, xIndex) in the
image interpolation */
public double precomputed_getWeightI(int l, int m, int u, int v)
{return prec_yWeight[v][l]*prec_xWeight[u][m];}
/*------------------------------------------------------------------*/
/* Interpolate the X and Y derivatives of the image at a
given point. */
public void precomputed_interpolateD(double []D, int u, int v) {
// Only SplineDegree=3 is implemented
D[0]=D[1]=0.0F;
for (int j = 0; j<4; j++) {
double sx=0.0F, sy=0.0F;
int iy=prec_yIndex[v][j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=prec_xIndex[u][i];
if (ix!=-1) {
double c;
if (fromCurrent) c=currentCoefficient[p + ix];
else c=coefficient[p + ix];
sx += prec_dxWeight[u][i]*c;
sy += prec_xWeight[u][i]*c;
}
}
D[0]+= prec_yWeight[v][j] * sx;
D[1]+=prec_dyWeight[v][j] * sy;
}
}
} /* end Interpolate D */
/*------------------------------------------------------------------*/
/* Interpolate the XY, XX and YY derivatives of the image at a
given point. */
public void precomputed_interpolateD2 (double []D2, int u, int v) {
// Only SplineDegree=3 is implemented
D2[0]=D2[1]=D2[2]=0.0F;
for (int j = 0; j<4; j++) {
double sxy=0.0F, sxx=0.0F, syy=0.0F;
int iy=prec_yIndex[v][j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=prec_xIndex[u][i];
if (ix!=-1) {
double c;
if (fromCurrent) c=currentCoefficient[p + ix];
else c=coefficient[p + ix];
sxy += prec_dxWeight[u][i]*c;
sxx += prec_d2xWeight[u][i]*c;
syy += prec_xWeight[u][i]*c;
}
}
D2[0]+= prec_dyWeight[v][j] * sxy;
D2[1]+= prec_yWeight[v][j] * sxx;
D2[2]+=prec_d2yWeight[v][j] * syy;
}
}
} /* end Interpolate dxdy, dxdx and dydy */
/*------------------------------------------------------------------*/
/* Interpolate the image at a given point. */
public double precomputed_interpolateI (int u, int v) {
// Only SplineDegree=3 is implemented
double ival=0.0F;
for (int j = 0; j<4; j++) {
double s=0.0F;
int iy=prec_yIndex[v][j];
if (iy!=-1) {
int p=iy*widthToUse;
for (int i=0; i<4; i++) {
int ix=prec_xIndex[u][i];
if (ix!=-1)
if (fromCurrent) s += prec_xWeight[u][i]*currentCoefficient[p + ix];
else s += prec_xWeight[u][i]*coefficient[p + ix];
}
ival+=prec_yWeight[v][j] * s;
}
}
return ival;
} /* end Interpolate Image */
/*------------------------------------------------------------------*/
/* Prepare precomputations for a given image size. */
public void precomputed_prepareForInterpolation(
int Ydim, int Xdim, int intervals) {
// Ask for memory
prec_xIndex =new int [Xdim][4];
prec_yIndex =new int [Ydim][4];
prec_xWeight =new double[Xdim][4];
prec_yWeight =new double[Ydim][4];
prec_dxWeight =new double[Xdim][4];
prec_dyWeight =new double[Ydim][4];
prec_d2xWeight=new double[Xdim][4];
prec_d2yWeight=new double[Ydim][4];
boolean ORIGINAL = false;
// Fill the precomputed weights and indexes for the Y axis
for (int v=0; v<Ydim; v++) {
// Express the current point in Spline units
final double tv = (double)(v * intervals) / (double)(Ydim - 1) + 1.0F;
final double tu = 1.0F;
// Compute all weights and indexes
prepareForInterpolation(tu,tv,ORIGINAL);
// Copy all values
for (int k=0; k<4; k++) {
prec_yIndex [v][k]= yIndex [k];
prec_yWeight [v][k]= yWeight[k];
prec_dyWeight [v][k]= dyWeight[k];
prec_d2yWeight[v][k]=d2yWeight[k];
}
}
// Fill the precomputed weights and indexes for the X axis
for (int u=0; u<Xdim; u++) {
// Express the current point in Spline units
final double tv = 1.0F;
final double tu = (double)(u * intervals) / (double)(Xdim - 1) + 1.0F;
// Compute all weights and indexes
prepareForInterpolation(tu,tv,ORIGINAL);
// Copy all values
for (int k=0; k<4; k++) {
prec_xIndex [u][k]= xIndex [k];
prec_xWeight [u][k]= xWeight[k];
prec_dxWeight [u][k]= dxWeight[k];
prec_d2xWeight[u][k]=d2xWeight[k];
}
}
}
/*------------------------------------------------------------------*/
/* Start the image precomputations. The computation of the B-spline
coefficients of the full-size image is not interruptible; all other
methods are. */
public void run (
) {
coefficient = getBasicFromCardinal2D();
buildCoefficientPyramid();
if (isTarget) buildImagePyramid();
} /* end run */
/*------------------------------------------------------------------*/
/* Set spline coefficients */
public void setCoefficients (
final double []c, // Set of B-spline coefficients
final int Ydim, // Dimensions of the set of coefficients
final int Xdim,
final int offset // Offset of the beginning of the array
// with respect to the origin of c
) {
// Copy the array of coefficients
System.arraycopy(c, offset, coefficient, 0, Ydim*Xdim);
}
/*------------------------------------------------------------------*/
/* Sets the depth up to which the pyramids should be computed. */
public void setPyramidDepth (
final int pyramidDepth
) {
int proposedPyramidDepth=pyramidDepth;
// Check what is the maximum depth allowed by the image
int currentWidth=width;
int currentHeight=height;
int scale=0;
while (currentWidth>=min_image_size && currentHeight>=min_image_size) {
currentWidth/=2;
currentHeight/=2;
scale++;
}
scale--;
if (proposedPyramidDepth>scale) proposedPyramidDepth=scale;
this.pyramidDepth = proposedPyramidDepth;
} /* end setPyramidDepth */
/*------------------------------------------------------------------*/
/* Converts the pixel array of the incoming ImageProcessor
object into a local double array. The flag is target enables the
computation of the derivative or not. */
public unwarpJImageModel (
final ImageProcessor ip,
final boolean isTarget
) {
// Initialize thread
t = new Thread(this);
t.setDaemon(true);
// Get image information
this.isTarget = isTarget;
width = ip.getWidth();
height = ip.getHeight();
twiceWidth = 2*width;
twiceHeight = 2*height;
coefficientsAreMirrored = true;
// Copy the pixel array
int k = 0;
image = new double[width * height];
unwarpJMiscTools.extractImage(ip,image);
// Resize the speedup arrays
xIndex = new int[4];
yIndex = new int[4];
xWeight = new double[4];
yWeight = new double[4];
dxWeight = new double[4];
dyWeight = new double[4];
d2xWeight = new double[4];
d2yWeight = new double[4];
} /* end unwarpJImage */
/* The same as before, but take the image from an array */
public unwarpJImageModel (
final double [][]img,
final boolean isTarget
) {
// Initialize thread
t = new Thread(this);
t.setDaemon(true);
// Get image information
this.isTarget = isTarget;
width = img[0].length;
height = img.length;
twiceWidth = 2*width;
twiceHeight = 2*height;
coefficientsAreMirrored = true;
// Copy the pixel array
int k = 0;
image = new double[width * height];
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
image[k] = img[y][x];
// Resize the speedup arrays
xIndex = new int[4];
yIndex = new int[4];
xWeight = new double[4];
yWeight = new double[4];
dxWeight = new double[4];
dyWeight = new double[4];
d2xWeight = new double[4];
d2yWeight = new double[4];
} /* end unwarpJImage */
/*------------------------------------------------------------------*/
/* Initialize the model from a set of coefficients */
public unwarpJImageModel (
final double [][]c // Set of B-spline coefficients
) {
// Get the size of the input array
currentHeight = height = c.length;
currentWidth = width = c[0].length;
twiceCurrentHeight = twiceHeight = 2*height;
twiceCurrentWidth = twiceWidth = 2*width;
coefficientsAreMirrored = false;
// Copy the array of coefficients
coefficient=new double[height*width];
int k=0;
for (int y=0; y<height; y++, k+= width)
System.arraycopy(c[y], 0, coefficient, k, width);
// Resize the speedup arrays
xIndex = new int[4];
yIndex = new int[4];
xWeight = new double[4];
yWeight = new double[4];
dxWeight = new double[4];
dyWeight = new double[4];
d2xWeight = new double[4];
d2yWeight = new double[4];
}
/*------------------------------------------------------------------*/
/* Initialize the model from a set of coefficients.
The same as the previous function but now the coefficients
are in a single row. */
public unwarpJImageModel (
final double []c, // Set of B-spline coefficients
final int Ydim, // Dimensions of the set of coefficients
final int Xdim,
final int offset // Offset of the beginning of the array
// with respect to the origin of c
) {
// Get the size of the input array
currentHeight = height = Ydim;
currentWidth = width = Xdim;
twiceCurrentHeight = twiceHeight = 2*height;
twiceCurrentWidth = twiceWidth = 2*width;
coefficientsAreMirrored = false;
// Copy the array of coefficients
coefficient=new double[height*width];
System.arraycopy(c, offset, coefficient, 0, height*width);
// Resize the speedup arrays
xIndex = new int[4];
yIndex = new int[4];
xWeight = new double[4];
yWeight = new double[4];
dxWeight = new double[4];
dyWeight = new double[4];
d2xWeight = new double[4];
d2yWeight = new double[4];
}
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void antiSymmetricFirMirrorOffBounds1D (
final double[] h,
final double[] c,
final double[] s
) {
if (2 <= c.length) {
s[0] = h[1] * (c[1] - c[0]);
for (int i = 1; (i < (s.length - 1)); i++) {
s[i] = h[1] * (c[i + 1] - c[i - 1]);
}
s[s.length - 1] = h[1] * (c[c.length - 1] - c[c.length - 2]);
} else s[0] = 0.0;
} /* end antiSymmetricFirMirrorOffBounds1D */
/*------------------------------------------------------------------*/
private void basicToCardinal2D (
final double[] basic,
final double[] cardinal,
final int width,
final int height,
final int degree
) {
final double[] hLine = new double[width];
final double[] vLine = new double[height];
final double[] hData = new double[width];
final double[] vData = new double[height];
double[] h = null;
switch (degree) {
case 3:
h = new double[2];
h[0] = 2.0 / 3.0;
h[1] = 1.0 / 6.0;
break;
case 7:
h = new double[4];
h[0] = 151.0 / 315.0;
h[1] = 397.0 / 1680.0;
h[2] = 1.0 / 42.0;
h[3] = 1.0 / 5040.0;
break;
default:
h = new double[1];
h[0] = 1.0;
}
for (int y = 0; ((y < height) && (!t.isInterrupted())); y++) {
extractRow(basic, y, hLine);
symmetricFirMirrorOffBounds1D(h, hLine, hData);
putRow(cardinal, y, hData);
}
for (int x = 0; ((x < width) && (!t.isInterrupted())); x++) {
extractColumn(cardinal, width, x, vLine);
symmetricFirMirrorOffBounds1D(h, vLine, vData);
putColumn(cardinal, width, x, vData);
}
} /* end basicToCardinal2D */
/*------------------------------------------------------------------*/
private void buildCoefficientPyramid (
) {
int fullWidth;
int fullHeight;
double[] fullDual = new double[width * height];
int halfWidth = width;
int halfHeight = height;
basicToCardinal2D(coefficient, fullDual, width, height, 7);
for (int depth = 1; ((depth <= pyramidDepth) && (!t.isInterrupted())); depth++) {
fullWidth = halfWidth;
fullHeight = halfHeight;
halfWidth /= 2;
halfHeight /= 2;
final double[] halfDual = getHalfDual2D(fullDual, fullWidth, fullHeight);
final double[] halfCoefficient = getBasicFromCardinal2D(halfDual, halfWidth, halfHeight, 7);
cpyramid.push(halfCoefficient);
cpyramid.push(new Integer(halfHeight));
cpyramid.push(new Integer(halfWidth));
fullDual = halfDual;
}
smallestWidth = halfWidth;
smallestHeight = halfHeight;
currentDepth=pyramidDepth+1;
} /* end buildCoefficientPyramid */
/*------------------------------------------------------------------*/
private void buildImagePyramid (
) {
int fullWidth;
int fullHeight;
double[] fullDual = new double[width * height];
int halfWidth = width;
int halfHeight = height;
cardinalToDual2D(image, fullDual, width, height, 3);
for (int depth = 1; ((depth <= pyramidDepth) && (!t.isInterrupted())); depth++) {
fullWidth = halfWidth;
fullHeight = halfHeight;
halfWidth /= 2;
halfHeight /= 2;
final double[] halfDual = getHalfDual2D(fullDual, fullWidth, fullHeight);
final double[] halfImage = new double[halfWidth * halfHeight];
dualToCardinal2D(halfDual, halfImage, halfWidth, halfHeight, 3);
imgpyramid.push(halfImage);
imgpyramid.push(new Integer(halfHeight));
imgpyramid.push(new Integer(halfWidth));
fullDual = halfDual;
}
} /* end buildImagePyramid */
/*------------------------------------------------------------------*/
private void cardinalToDual2D (
final double[] cardinal,
final double[] dual,
final int width,
final int height,
final int degree
) {
basicToCardinal2D(getBasicFromCardinal2D(cardinal, width, height, degree),
dual, width, height, 2 * degree + 1);
} /* end cardinalToDual2D */
/*------------------------------------------------------------------*/
private void coefficientToGradient1D (
final double[] c
) {
final double[] h = {0.0, 1.0 / 2.0};
final double[] s = new double[c.length];
antiSymmetricFirMirrorOffBounds1D(h, c, s);
System.arraycopy(s, 0, c, 0, s.length);
} /* end coefficientToGradient1D */
/*------------------------------------------------------------------*/
private void coefficientToSamples1D (
final double[] c
) {
final double[] h = {2.0 / 3.0, 1.0 / 6.0};
final double[] s = new double[c.length];
symmetricFirMirrorOffBounds1D(h, c, s);
System.arraycopy(s, 0, c, 0, s.length);
} /* end coefficientToSamples1D */
/*------------------------------------------------------------------*/
private void coefficientToXYGradient2D (
final double[] basic,
final double[] xGradient,
final double[] yGradient,
final int width,
final int height
) {
final double[] hLine = new double[width];
final double[] hData = new double[width];
final double[] vLine = new double[height];
for (int y = 0; ((y < height) && (!t.isInterrupted())); y++) {
extractRow(basic, y, hLine);
System.arraycopy(hLine, 0, hData, 0, width);
coefficientToGradient1D(hLine);
coefficientToSamples1D(hData);
putRow(xGradient, y, hLine);
putRow(yGradient, y, hData);
}
for (int x = 0; ((x < width) && (!t.isInterrupted())); x++) {
extractColumn(xGradient, width, x, vLine);
coefficientToSamples1D(vLine);
putColumn(xGradient, width, x, vLine);
extractColumn(yGradient, width, x, vLine);
coefficientToGradient1D(vLine);
putColumn(yGradient, width, x, vLine);
}
} /* end coefficientToXYGradient2D */
/*------------------------------------------------------------------*/
private void dualToCardinal2D (
final double[] dual,
final double[] cardinal,
final int width,
final int height,
final int degree
) {
basicToCardinal2D(getBasicFromCardinal2D(dual, width, height, 2 * degree + 1),
cardinal, width, height, degree);
} /* end dualToCardinal2D */
/*------------------------------------------------------------------*/
private void extractColumn (
final double[] array,
final int width,
int x,
final double[] column
) {
for (int i = 0; (i < column.length); i++, x+=width)
column[i] = (double)array[x];
} /* end extractColumn */
/*------------------------------------------------------------------*/
private void extractRow (
final double[] array,
int y,
final double[] row
) {
y *= row.length;
for (int i = 0; (i < row.length); i++)
row[i] = (double)array[y++];
} /* end extractRow */
/*------------------------------------------------------------------*/
private double[] getBasicFromCardinal2D (
) {
final double[] basic = new double[width * height];
final double[] hLine = new double[width];
final double[] vLine = new double[height];
for (int y = 0; (y < height); y++) {
extractRow(image, y, hLine);
samplesToInterpolationCoefficient1D(hLine, 3, 0.0);
putRow(basic, y, hLine);
}
for (int x = 0; (x < width); x++) {
extractColumn(basic, width, x, vLine);
samplesToInterpolationCoefficient1D(vLine, 3, 0.0);
putColumn(basic, width, x, vLine);
}
return(basic);
} /* end getBasicFromCardinal2D */
/*------------------------------------------------------------------*/
private double[] getBasicFromCardinal2D (
final double[] cardinal,
final int width,
final int height,
final int degree
) {
final double[] basic = new double[width * height];
final double[] hLine = new double[width];
final double[] vLine = new double[height];
for (int y = 0; ((y < height) && (!t.isInterrupted())); y++) {
extractRow(cardinal, y, hLine);
samplesToInterpolationCoefficient1D(hLine, degree, 0.0);
putRow(basic, y, hLine);
}
for (int x = 0; ((x < width) && (!t.isInterrupted())); x++) {
extractColumn(basic, width, x, vLine);
samplesToInterpolationCoefficient1D(vLine, degree, 0.0);
putColumn(basic, width, x, vLine);
}
return(basic);
} /* end getBasicFromCardinal2D */
/*------------------------------------------------------------------*/
private double[] getHalfDual2D (
final double[] fullDual,
final int fullWidth,
final int fullHeight
) {
final int halfWidth = fullWidth / 2;
final int halfHeight = fullHeight / 2;
final double[] hLine = new double[fullWidth];
final double[] hData = new double[halfWidth];
final double[] vLine = new double[fullHeight];
final double[] vData = new double[halfHeight];
final double[] demiDual = new double[halfWidth * fullHeight];
final double[] halfDual = new double[halfWidth * halfHeight];
for (int y = 0; ((y < fullHeight) && (!t.isInterrupted())); y++) {
extractRow(fullDual, y, hLine);
reduceDual1D(hLine, hData);
putRow(demiDual, y, hData);
}
for (int x = 0; ((x < halfWidth) && (!t.isInterrupted())); x++) {
extractColumn(demiDual, halfWidth, x, vLine);
reduceDual1D(vLine, vData);
putColumn(halfDual, halfWidth, x, vData);
}
return(halfDual);
} /* end getHalfDual2D */
/*------------------------------------------------------------------*/
private double getInitialAntiCausalCoefficientMirrorOffBounds (
final double[] c,
final double z,
final double tolerance
) {
return(z * c[c.length - 1] / (z - 1.0));
} /* end getInitialAntiCausalCoefficientMirrorOffBounds */
/*------------------------------------------------------------------*/
private double getInitialCausalCoefficientMirrorOffBounds (
final double[] c,
final double z,
final double tolerance
) {
double z1 = z;
double zn = Math.pow(z, c.length);
double sum = (1.0 + z) * (c[0] + zn * c[c.length - 1]);
int horizon = c.length;
if (0.0 < tolerance) {
horizon = 2 + (int)(Math.log(tolerance) / Math.log(Math.abs(z)));
horizon = (horizon < c.length) ? (horizon) : (c.length);
}
zn = zn * zn;
for (int n = 1; (n < (horizon - 1)); n++) {
z1 = z1 * z;
zn = zn / z;
sum = sum + (z1 + zn) * c[n];
}
return(sum / (1.0 - Math.pow(z, 2 * c.length)));
} /* end getInitialCausalCoefficientMirrorOffBounds */
/*------------------------------------------------------------------*/
private void putColumn (
final double[] array,
final int width,
int x,
final double[] column
) {
for (int i = 0; (i < column.length); i++, x+=width)
array[x] = (double)column[i];
} /* end putColumn */
/*------------------------------------------------------------------*/
private void putRow (
final double[] array,
int y,
final double[] row
) {
y *= row.length;
for (int i = 0; (i < row.length); i++)
array[y++] = (double)row[i];
} /* end putRow */
/*------------------------------------------------------------------*/
private void reduceDual1D (
final double[] c,
final double[] s
) {
final double h[] = {6.0 / 16.0, 4.0 / 16.0, 1.0 / 16.0};
if (2 <= s.length) {
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2]);
for (int i = 2, j = 1; (j < (s.length - 1)); i += 2, j++) {
s[j] = h[0] * c[i] + h[1] * (c[i - 1] + c[i + 1])
+ h[2] * (c[i - 2] + c[i + 2]);
}
if (c.length == (2 * s.length)) {
s[s.length - 1] = h[0] * c[c.length - 2] + h[1] * (c[c.length - 3] + c[c.length - 1])
+ h[2] * (c[c.length - 4] + c[c.length - 1]);
}
else {
s[s.length - 1] = h[0] * c[c.length - 3] + h[1] * (c[c.length - 4] + c[c.length - 2])
+ h[2] * (c[c.length - 5] + c[c.length - 1]);
}
}
else {
switch (c.length) {
case 3:
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2]);
break;
case 2:
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + 2.0 * h[2] * c[1];
break;
default:
}
}
} /* end reduceDual1D */
/*------------------------------------------------------------------*/
private void samplesToInterpolationCoefficient1D (
final double[] c,
final int degree,
final double tolerance
) {
double[] z = new double[0];
double lambda = 1.0;
switch (degree) {
case 3:
z = new double[1];
z[0] = Math.sqrt(3.0) - 2.0;
break;
case 7:
z = new double[3];
z[0] = -0.5352804307964381655424037816816460718339231523426924148812;
z[1] = -0.122554615192326690515272264359357343605486549427295558490763;
z[2] = -0.0091486948096082769285930216516478534156925639545994482648003;
break;
default:
}
if (c.length == 1) {
return;
}
for (int k = 0; (k < z.length); k++) {
lambda *= (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
}
for (int n = 0; (n < c.length); n++) {
c[n] = c[n] * lambda;
}
for (int k = 0; (k < z.length); k++) {
c[0] = getInitialCausalCoefficientMirrorOffBounds(c, z[k], tolerance);
for (int n = 1; (n < c.length); n++) {
c[n] = c[n] + z[k] * c[n - 1];
}
c[c.length - 1] = getInitialAntiCausalCoefficientMirrorOffBounds(c, z[k], tolerance);
for (int n = c.length - 2; (0 <= n); n--) {
c[n] = z[k] * (c[n+1] - c[n]);
}
}
} /* end samplesToInterpolationCoefficient1D */
/*------------------------------------------------------------------*/
private void symmetricFirMirrorOffBounds1D (
final double[] h,
final double[] c,
final double[] s
) {
switch (h.length) {
case 2:
if (2 <= c.length) {
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]);
for (int i = 1; (i < (s.length - 1)); i++) {
s[i] = h[0] * c[i] + h[1] * (c[i - 1] + c[i + 1]);
}
s[s.length - 1] = h[0] * c[c.length - 1]
+ h[1] * (c[c.length - 2] + c[c.length - 1]);
}
else {
s[0] = (h[0] + 2.0 * h[1]) * c[0];
}
break;
case 4:
if (6 <= c.length) {
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2])
+ h[3] * (c[2] + c[3]);
s[1] = h[0] * c[1] + h[1] * (c[0] + c[2]) + h[2] * (c[0] + c[3])
+ h[3] * (c[1] + c[4]);
s[2] = h[0] * c[2] + h[1] * (c[1] + c[3]) + h[2] * (c[0] + c[4])
+ h[3] * (c[0] + c[5]);
for (int i = 3; (i < (s.length - 3)); i++) {
s[i] = h[0] * c[i] + h[1] * (c[i - 1] + c[i + 1])
+ h[2] * (c[i - 2] + c[i + 2]) + h[3] * (c[i - 3] + c[i + 3]);
}
s[s.length - 3] = h[0] * c[c.length - 3]
+ h[1] * (c[c.length - 4] + c[c.length - 2])
+ h[2] * (c[c.length - 5] + c[c.length - 1])
+ h[3] * (c[c.length - 6] + c[c.length - 1]);
s[s.length - 2] = h[0] * c[c.length - 2]
+ h[1] * (c[c.length - 3] + c[c.length - 1])
+ h[2] * (c[c.length - 4] + c[c.length - 1])
+ h[3] * (c[c.length - 5] + c[c.length - 2]);
s[s.length - 1] = h[0] * c[c.length - 1]
+ h[1] * (c[c.length - 2] + c[c.length - 1])
+ h[2] * (c[c.length - 3] + c[c.length - 2])
+ h[3] * (c[c.length - 4] + c[c.length - 3]);
}
else {
switch (c.length) {
case 5:
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2])
+ h[3] * (c[2] + c[3]);
s[1] = h[0] * c[1] + h[1] * (c[0] + c[2]) + h[2] * (c[0] + c[3])
+ h[3] * (c[1] + c[4]);
s[2] = h[0] * c[2] + h[1] * (c[1] + c[3])
+ (h[2] + h[3]) * (c[0] + c[4]);
s[3] = h[0] * c[3] + h[1] * (c[2] + c[4]) + h[2] * (c[1] + c[4])
+ h[3] * (c[0] + c[3]);
s[4] = h[0] * c[4] + h[1] * (c[3] + c[4]) + h[2] * (c[2] + c[3])
+ h[3] * (c[1] + c[2]);
break;
case 4:
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2])
+ h[3] * (c[2] + c[3]);
s[1] = h[0] * c[1] + h[1] * (c[0] + c[2]) + h[2] * (c[0] + c[3])
+ h[3] * (c[1] + c[3]);
s[2] = h[0] * c[2] + h[1] * (c[1] + c[3]) + h[2] * (c[0] + c[3])
+ h[3] * (c[0] + c[2]);
s[3] = h[0] * c[3] + h[1] * (c[2] + c[3]) + h[2] * (c[1] + c[2])
+ h[3] * (c[0] + c[1]);
break;
case 3:
s[0] = h[0] * c[0] + h[1] * (c[0] + c[1]) + h[2] * (c[1] + c[2])
+ 2.0 * h[3] * c[2];
s[1] = h[0] * c[1] + (h[1] + h[2]) * (c[0] + c[2])
+ 2.0 * h[3] * c[1];
s[2] = h[0] * c[2] + h[1] * (c[1] + c[2]) + h[2] * (c[0] + c[1])
+ 2.0 * h[3] * c[0];
break;
case 2:
s[0] = (h[0] + h[1] + h[3]) * c[0] + (h[1] + 2.0 * h[2] + h[3]) * c[1];
s[1] = (h[0] + h[1] + h[3]) * c[1] + (h[1] + 2.0 * h[2] + h[3]) * c[0];
break;
case 1:
s[0] = (h[0] + 2.0 * (h[1] + h[2] + h[3])) * c[0];
break;
default:
}
}
break;
default:
}
} /* end symmetricFirMirrorOffBounds1D */
} /* end class unwarpJImageModel */
/*====================================================================
| unwarpJMask
\===================================================================*/
/* This class is responsible for the mask preprocessing that takes
place concurrently with user-interface events. It contains methods
to compute the mask pyramids. */
class unwarpJMask
{ /* begin class unwarpJMask */
/*....................................................................
Private variables
....................................................................*/
// Mask related
private boolean[] mask;
private int width;
private int height;
private Polygon polygon=null;
private boolean mask_from_the_stack;
/*....................................................................
Public methods
....................................................................*/
/* Bounding box for the mask.
An array is returned with the convention [x0,y0,xF,yF]. This array
is returned in corners. This vector should be already resized. */
public void BoundingBox(int [] corners) {
if (polygon.npoints!=0) {
Rectangle boundingbox=polygon.getBounds();
corners[0]=(int)boundingbox.x;
corners[1]=(int)boundingbox.y;
corners[2]=corners[0]+(int)boundingbox.width;
corners[3]=corners[1]+(int)boundingbox.height;
} else {
corners[0]=0;
corners[1]=0;
corners[2]=width;
corners[3]=height;
}
}
/*------------------------------------------------------------------*/
/* Set to true every pixel of the full-size mask. */
public void clearMask (
) {
int k = 0;
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++)
mask[k++] = true;
polygon=new Polygon();
} /* end clearMask */
/*------------------------------------------------------------------*/
/* Fill the mask associated to the mask points. */
public void fillMask (int tool) {
int k=0;
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++) {
mask[k] = polygon.contains(x,y);
if (tool==unwarpJPointAction.INVERTMASK) mask[k]=!mask[k];
k++;
}
}
/*------------------------------------------------------------------*/
/* Returns the value of the mask at a certain pixel.
If the sample is not integer then the closest point is returned. */
public boolean getValue(double x, double y) {
int u=(int)Math.round(x);
int v=(int)Math.round(y);
if (u<0 || u>=width || v<0 || v>=height) return false;
else return mask[v*width+u];
}
/*------------------------------------------------------------------*/
/* Get a point from the mask. */
public Point getPoint(int i) {
return new Point(polygon.xpoints[i],polygon.ypoints[i]);
}
/*------------------------------------------------------------------*/
/* True if the mask was taken from the stack. */
public boolean isFromStack() {
return mask_from_the_stack;
}
/*------------------------------------------------------------------*/
/* Get the number of points in the mask. */
public int numberOfMaskPoints() {return polygon.npoints;}
/*------------------------------------------------------------------*/
/* Read mask from file.
An error is shown if the file read is not of the same size as the
previous mask. */
public void readFile(String filename) {
ImagePlus aux = new ImagePlus(filename);
if (aux.getWidth()!=width || aux.getHeight()!=height)
IJ.error("Mask in file is not of the expected size");
ImageProcessor ip=aux.getProcessor();
int k=0;
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
if (ip.getPixelValue(x,y)!=0) mask[k]=true;
else mask[k]=false;
}
/*------------------------------------------------------------------*/
/* Show mask. */
public void showMask () {
double [][]img=new double[height][width];
int k = 0;
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++)
if (mask[k++]) img[y][x]=1; else img[y][x]=0;
unwarpJMiscTools.showImage("Mask",img);
}
/*------------------------------------------------------------------*/
/* Set the mask points. */
public void setMaskPoints (final Vector listMaskPoints) {
int imax=listMaskPoints.size();
for (int i=0; i<imax; i++) {
Point p=(Point)listMaskPoints.elementAt(i);
polygon.addPoint(p.x,p.y);
}
}
/*------------------------------------------------------------------*/
/* Sets the value of the mask at a certain pixel. */
public void setValue(int u, int v, boolean value) {
if (u>=0 && u<width && v>=0 && v<height) mask[v*width+u]=value;
}
/*------------------------------------------------------------------*/
/* Empty constructor, the input image is used only to take the
image size. */
public unwarpJMask (
final ImageProcessor ip, boolean take_mask
) {
width = ip.getWidth();
height = ip.getHeight();
mask = new boolean[width * height];
if (!take_mask) {
mask_from_the_stack=false;
clearMask();
} else {
mask_from_the_stack=true;
int k=0;
if (ip instanceof ByteProcessor) {
final byte[] pixels = (byte[])ip.getPixels();
for (int y = 0; (y < height); y++) {
for (int x = 0; (x < width); x++, k++) {
mask[k] = (pixels[k] != 0);
}
}
}
else if (ip instanceof ShortProcessor) {
final short[] pixels = (short[])ip.getPixels();
for (int y = 0; (y < height); y++) {
for (int x = 0; (x < width); x++, k++) {
mask[k] = (pixels[k] != 0);
}
}
}
else if (ip instanceof FloatProcessor) {
final float[] pixels = (float[])ip.getPixels();
for (int y = 0; (y < height); y++) {
for (int x = 0; (x < width); x++, k++) {
mask[k] = (pixels[k] != 0.0F);
}
}
}
}
} /* end unwarpJMask */
} /* end class unwarpJMask */
/*====================================================================
| unwarpJMathTools
\===================================================================*/
class unwarpJMathTools {
private static final double FLT_EPSILON = (double)Float.intBitsToFloat((int)0x33FFFFFF);
private static final int MAX_SVD_ITERATIONS = 1000;
/*------------------------------------------------------------------*/
public static double Bspline01 (
double x
) {
x = Math.abs(x);
if (x < 1.0F) {
return(1.0F - x);
}
else {
return(0.0F);
}
} /* Bspline01 */
/*------------------------------------------------------------------*/
public static double Bspline02 (
double x
) {
x = Math.abs(x);
if (x < 0.5F) {
return(3.0F / 4.0F - x * x);
}
else if (x < 1.5F) {
x -= 1.5F;
return(0.5F * x * x);
}
else {
return(0.0F);
}
} /* Bspline02 */
/*------------------------------------------------------------------*/
public static double Bspline03 (
double x
) {
x = Math.abs(x);
if (x < 1.0F) {
return(0.5F * x * x * (x - 2.0F) + (2.0F / 3.0F));
}
else if (x < 2.0F) {
x -= 2.0F;
return(x * x * x / (-6.0F));
}
else {
return(0.0F);
}
} /* Bspline03 */
/*------------------------------------------------------------------*/
public static double EuclideanNorm (
final double a,
final double b
) {
final double absa = Math.abs(a);
final double absb = Math.abs(b);
if (absb < absa) {
return(absa * Math.sqrt(1.0 + (absb * absb / (absa * absa))));
}
else {
return((absb == 0.0F) ? (0.0F)
: (absb * Math.sqrt(1.0 + (absa * absa / (absb * absb)))));
}
} /* end EuclideanNorm */
/*------------------------------------------------------------------*/
public static boolean invertMatrixSVD(
int Ydim, // Input,
int Xdim, // Input,
final double [][]B, // Input, matrix to invert
final double [][]iB // Output, inverted matrix
) {
boolean underconstrained=false;
final double[] W = new double[Xdim];
final double[][] V = new double[Xdim][Xdim];
// B=UWV^t (U is stored in B)
singularValueDecomposition(B, W, V);
// B^-1=VW^-1U^t
// Compute W^-1
int Nzeros=0;
for (int k = 0; k<Xdim; k++) {
if (Math.abs(W[k]) < FLT_EPSILON) {
W[k] = 0.0F;
Nzeros++;
} else
W[k] = 1.0F / W[k];
}
if (Ydim-Nzeros<Xdim) underconstrained=true;
// Compute VW^-1
for (int i = 0; i<Xdim; i++)
for (int j = 0; j<Xdim; j++)
V[i][j] *= W[j];
// Compute B^-1
// iB should have been already resized
for (int i = 0; i<Xdim; i++) {
for (int j = 0; j<Ydim; j++) {
iB[i][j] = 0.0F;
for (int k = 0; k<Xdim; k++)
iB[i][j] += V[i][k] * B[j][k];
}
}
return underconstrained;
} /* invertMatrixSVD */
/*------------------------------------------------------------------*/
/*********************************************************************
* Gives the least-squares solution to (A * x = b) such that
* (A^T * A)^-1 * A^T * b = x is a vector of size (column), where A is
* a (line x column) matrix, and where b is a vector of size (line).
* The result may differ from that obtained by a singular-value
* decomposition in the cases where the least-squares solution is not
* uniquely defined (SVD returns the solution of least norm, not QR).
*
* @param A An input matrix A[line][column] of size (line x column)
* @param b An input vector b[line] of size (line)
* @return An output vector x[column] of size (column)
********************************************************************/
public static 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.0F;
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.0F) {
x[i] = 0.0F;
}
else {
x[i] /= s;
}
for (int j = i - 1; (0 <= j); j--) {
x[j] -= R[j][i] * x[i];
}
}
return(x);
} /* end linearLeastSquares */
/*------------------------------------------------------------------*/
public static double nchoosek(int n, int k) {
if (k>n) return 0;
if (k==0) return 1;
if (k==1) return n;
if (k>n/2) k=n-k;
double prod=1;
for (int i=1; i<=k; i++) prod*=(n-k+i)/i;
return prod;
}
/*------------------------------------------------------------------*/
/*********************************************************************
* Decomposes the (line x column) input matrix Q into an orthonormal
* output matrix Q of same size (line x column) and an upper-diagonal
* square matrix R of size (column x column), such that the matrix
* product (Q * R) gives the input matrix, and such that the matrix
* product (Q^T * Q) gives the identity.
*
* @param Q An in-place (line x column) matrix Q[line][column], which
* expects as input the matrix to decompose, and which returns as
* output an orthonormal matrix
* @param R An output (column x column) square matrix R[column][column]
********************************************************************/
public static 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.0F;
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.0F;
for (int i = 0; (i < lines); i++) {
s += Q[i][j] * Q[i][j];
}
if ((s * s) == 0.0F) {
s = 0.0F;
}
else {
s = 1.0F / 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.0F;
}
for (int j = i; (j < columns); j++) {
R[i][j] = 0.0F;
for (int k = 0; (k < lines); k++) {
R[i][j] += Q[k][i] * A[k][j];
}
}
}
} /* end QRdecomposition */
/* -----------------------------------------------------------------*/
public static void showMatrix(
int Ydim,
int Xdim,
final double [][]A
) {
for (int i=0; i<Ydim; i++) {
for (int j=0; j<Xdim; j++)
System.out.print(A[i][j]+" ");
System.out.println();
}
}
/*------------------------------------------------------------------*/
public static void singularValueDecomposition (
final double[][] U,
final double[] W,
final double[][] V
) {
final int lines = U.length;
final int columns = U[0].length;
final double[] rv1 = new double[columns];
double norm, scale;
double c, f, g, h, s;
double x, y, z;
int l = 0;
int nm = 0;
boolean flag;
g = scale = norm = 0.0F;
for (int i = 0; (i < columns); i++) {
l = i + 1;
rv1[i] = scale * g;
g = s = scale = 0.0F;
if (i < lines) {
for (int k = i; (k < lines); k++) {
scale += Math.abs(U[k][i]);
}
if (scale != 0.0) {
for (int k = i; (k < lines); k++) {
U[k][i] /= scale;
s += U[k][i] * U[k][i];
}
f = U[i][i];
g = (0.0 <= f) ? (-Math.sqrt((double)s))
: (Math.sqrt((double)s));
h = f * g - s;
U[i][i] = f - g;
for (int j = l; (j < columns); j++) {
s = 0.0F;
for (int k = i; (k < lines); k++) {
s += U[k][i] * U[k][j];
}
f = s / h;
for (int k = i; (k < lines); k++) {
U[k][j] += f * U[k][i];
}
}
for (int k = i; (k < lines); k++) {
U[k][i] *= scale;
}
}
}
W[i] = scale * g;
g = s = scale = 0.0F;
if ((i < lines) && (i != (columns - 1))) {
for (int k = l; (k < columns); k++) {
scale += Math.abs(U[i][k]);
}
if (scale != 0.0) {
for (int k = l; (k < columns); k++) {
U[i][k] /= scale;
s += U[i][k] * U[i][k];
}
f = U[i][l];
g = (0.0 <= f) ? (-Math.sqrt(s))
: (Math.sqrt(s));
h = f * g - s;
U[i][l] = f - g;
for (int k = l; (k < columns); k++) {
rv1[k] = U[i][k] / h;
}
for (int j = l; (j < lines); j++) {
s = 0.0F;
for (int k = l; (k < columns); k++) {
s += U[j][k] * U[i][k];
}
for (int k = l; (k < columns); k++) {
U[j][k] += s * rv1[k];
}
}
for (int k = l; (k < columns); k++) {
U[i][k] *= scale;
}
}
}
norm = ((Math.abs(W[i]) + Math.abs(rv1[i])) < norm) ? (norm)
: (Math.abs(W[i]) + Math.abs(rv1[i]));
}
for (int i = columns - 1; (0 <= i); i--) {
if (i < (columns - 1)) {
if (g != 0.0) {
for (int j = l; (j < columns); j++) {
V[j][i] = U[i][j] / (U[i][l] * g);
}
for (int j = l; (j < columns); j++) {
s = 0.0F;
for (int k = l; (k < columns); k++) {
s += U[i][k] * V[k][j];
}
for (int k = l; (k < columns); k++) {
if (s != 0.0) {
V[k][j] += s * V[k][i];
}
}
}
}
for (int j = l; (j < columns); j++) {
V[i][j] = V[j][i] = 0.0F;
}
}
V[i][i] = 1.0F;
g = rv1[i];
l = i;
}
for (int i = (lines < columns) ? (lines - 1) : (columns - 1); (0 <= i); i--) {
l = i + 1;
g = W[i];
for (int j = l; (j < columns); j++) {
U[i][j] = 0.0F;
}
if (g != 0.0) {
g = 1.0F / g;
for (int j = l; (j < columns); j++) {
s = 0.0F;
for (int k = l; (k < lines); k++) {
s += U[k][i] * U[k][j];
}
f = s * g / U[i][i];
for (int k = i; (k < lines); k++) {
if (f != 0.0) {
U[k][j] += f * U[k][i];
}
}
}
for (int j = i; (j < lines); j++) {
U[j][i] *= g;
}
}
else {
for (int j = i; (j < lines); j++) {
U[j][i] = 0.0F;
}
}
U[i][i] += 1.0F;
}
for (int k = columns - 1; (0 <= k); k--) {
for (int its = 1; (its <= MAX_SVD_ITERATIONS); its++) {
flag = true;
for (l = k; (0 <= l); l--) {
nm = l - 1;
if ((Math.abs(rv1[l]) + norm) == norm) {
flag = false;
break;
}
if ((Math.abs(W[nm]) + norm) == norm) {
break;
}
}
if (flag) {
c = 0.0F;
s = 1.0F;
for (int i = l; (i <= k); i++) {
f = s * rv1[i];
rv1[i] *= c;
if ((Math.abs(f) + norm) == norm) {
break;
}
g = W[i];
h = EuclideanNorm(f, g);
W[i] = h;
h = 1.0F / h;
c = g * h;
s = -f * h;
for (int j = 0; (j < lines); j++) {
y = U[j][nm];
z = U[j][i];
U[j][nm] = y * c + z * s;
U[j][i] = z * c - y * s;
}
}
}
z = W[k];
if (l == k) {
if (z < 0.0) {
W[k] = -z;
for (int j = 0; (j < columns); j++) {
V[j][k] = -V[j][k];
}
}
break;
}
if (its == MAX_SVD_ITERATIONS) {
return;
}
x = W[l];
nm = k - 1;
y = W[nm];
g = rv1[nm];
h = rv1[k];
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0F * h * y);
g = EuclideanNorm(f, 1.0F);
f = ((x - z) * (x + z) + h * ((y / (f + ((0.0 <= f) ? (Math.abs(g))
: (-Math.abs(g))))) - h)) / x;
c = s = 1.0F;
for (int j = l; (j <= nm); j++) {
int i = j + 1;
g = rv1[i];
y = W[i];
h = s * g;
g = c * g;
z = EuclideanNorm(f, h);
rv1[j] = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = g * c - x * s;
h = y * s;
y *= c;
for (int jj = 0; (jj < columns); jj++) {
x = V[jj][j];
z = V[jj][i];
V[jj][j] = x * c + z * s;
V[jj][i] = z * c - x * s;
}
z = EuclideanNorm(f, h);
W[j] = z;
if (z != 0.0F) {
z = 1.0F / z;
c = f * z;
s = h * z;
}
f = c * g + s * y;
x = c * y - s * g;
for (int jj = 0; (jj < lines); jj++) {
y = U[jj][j];
z = U[jj][i];
U[jj][j] = y * c + z * s;
U[jj][i] = z * c - y * s;
}
}
rv1[l] = 0.0F;
rv1[k] = f;
W[k] = x;
}
}
} /* end singularValueDecomposition */
public static void singularValueBackSubstitution (
final double [][]U, /* input matrix */
final double []W, /* vector of singular values */
final double [][]V, /* untransposed orthogonal matrix */
final double []B, /* input vector */
final double []X /* returned solution */
){
/* solve (U.W.Transpose(V)).X == B in terms of X */
/* {U, W, V} are given by SingularValueDecomposition */
/* by convention, set w[i,j]=0 to get (1/w[i,j])=0 */
/* the size of the input matrix U is (Lines x Columns) */
/* the size of the vector (1/W) of singular values is (Columns) */
/* the size of the untransposed orthogonal matrix V is (Columns x Columns) */
/* the size of the input vector B is (Lines) */
/* the size of the output vector X is (Columns) */
final int Lines = U.length;
final int Columns = U[0].length;
double [] aux = new double [Columns];
// A=UWV^t
// A^-1=VW^-1U^t
// X=A^-1*B
// Perform aux=W^-1 U^t B
for (int i=0; i<Columns; i++) {
aux[i]=0.0F;
if (Math.abs(W[i]) > FLT_EPSILON) {
for (int j=0; j<Lines; j++) aux[i]+=U[j][i]*B[j]; // U^t B
aux[i]/=W[i]; // W^-1 U^t B
}
}
// Perform X=V aux
for (int i=0; i<Columns; i++) {
X[i]=0.0F;
for (int j=0; j<Columns; j++) X[i]+=V[i][j]*aux[j];
}
}
} /* End MathTools */
/*====================================================================
| unwarpJMiscTools
\===================================================================*/
class unwarpJMiscTools {
/* Apply a given splines transformation to the source image.
The source image is modified. The target image is used to know
the output size. */
static public void applyTransformationToSource(
ImagePlus sourceImp, ImagePlus targetImp,
unwarpJImageModel source,
int intervals,
double [][]cx,
double [][]cy) {
int targetHeight=targetImp.getProcessor().getHeight();
int targetWidth =targetImp.getProcessor().getWidth ();
int sourceHeight=sourceImp.getProcessor().getHeight();
int sourceWidth =sourceImp.getProcessor().getWidth ();
// Ask for memory for the transformation
double [][] transformation_x=new double [targetHeight][targetWidth];
double [][] transformation_y=new double [targetHeight][targetWidth];
// Compute the deformation
// Set these coefficients to an interpolator
unwarpJImageModel swx = new unwarpJImageModel(cx);
unwarpJImageModel swy = new unwarpJImageModel(cy);
// Compute the transformation mapping
boolean ORIGINAL=false;
for (int v=0; v<targetHeight; v++) {
final double tv = (double)(v * intervals) / (double)(targetHeight - 1) + 1.0F;
for (int u = 0; u<targetWidth; u++) {
final double tu = (double)(u * intervals) / (double)(targetWidth - 1) + 1.0F;
swx.prepareForInterpolation(tu, tv, ORIGINAL); transformation_x[v][u] = swx.interpolateI();
swy.prepareForInterpolation(tu, tv, ORIGINAL); transformation_y[v][u] = swy.interpolateI();
}
}
// Compute the warped image
FloatProcessor fp=new FloatProcessor(targetWidth,targetHeight);
for (int v=0; v<targetHeight; v++)
for (int u=0; u<targetWidth; u++) {
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
if (x>=0 && x<sourceWidth && y>=0 && y<sourceHeight) {
source.prepareForInterpolation(x,y,ORIGINAL);
fp.putPixelValue(u,v,source.interpolateI());
} else
fp.putPixelValue(u,v,0);
}
fp.resetMinAndMax();
sourceImp.setProcessor(sourceImp.getTitle(),fp);
sourceImp.updateImage();
}
/*------------------------------------------------------------------*/
/* Draw an arrow between two points.
The arrow head is in (x2,y2) */
static public void drawArrow(double [][]canvas, int x1, int y1,
int x2, int y2, double color, int arrow_size) {
drawLine(canvas,x1,y1,x2,y2,color);
int arrow_size2=2*arrow_size;
// Do not draw the arrow_head if the arrow is very small
if ((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)<arrow_size*arrow_size) return;
// Vertical arrow
if (x2 == x1) {
if (y2 > y1) {
drawLine(canvas,x2,y2,x2-arrow_size,y2-arrow_size2,color);
drawLine(canvas,x2,y2,x2+arrow_size,y2-arrow_size2,color);
} else {
drawLine(canvas,x2,y2,x2-arrow_size,y2+arrow_size2,color);
drawLine(canvas,x2,y2,x2+arrow_size,y2+arrow_size2,color);
}
}
// Horizontal arrow
else if (y2 == y1) {
if (x2 > x1) {
drawLine(canvas,x2,y2,x2-arrow_size2,y2-arrow_size,color);
drawLine(canvas,x2,y2,x2-arrow_size2,y2+arrow_size,color);
} else {
drawLine(canvas,x2,y2,x2+arrow_size2,y2-arrow_size,color);
drawLine(canvas,x2,y2,x2+arrow_size2,y2+arrow_size,color);
}
}
// Now we need to rotate the arrow head about the origin
else {
// Calculate the angle of rotation and adjust for the quadrant
double t1 = Math.abs(new Integer(y2 - y1).doubleValue());
double t2 = Math.abs(new Integer(x2 - x1).doubleValue());
double theta = Math.atan(t1 / t2);
if (x2 < x1) {
if (y2 < y1) theta = Math.PI + theta;
else theta = - (Math.PI + theta);
} else if (x2 > x1 && y2 < y1)
theta = 2*Math.PI - theta;
double cosTheta = Math.cos(theta);
double sinTheta = Math.sin(theta);
// Create the other points and translate the arrow to the origin
Point p2 = new Point(-arrow_size2,-arrow_size);
Point p3 = new Point(-arrow_size2,+arrow_size);
// Rotate the points (without using matrices!)
int x = new Long(Math.round((cosTheta * p2.x) - (sinTheta * p2.y))).intValue();
p2.y = new Long(Math.round((sinTheta * p2.x) + (cosTheta * p2.y))).intValue();
p2.x = x;
x = new Long(Math.round((cosTheta * p3.x) - (sinTheta * p3.y))).intValue();
p3.y = new Long(Math.round((sinTheta * p3.x) + (cosTheta * p3.y))).intValue();
p3.x = x;
// Translate back to desired location and add to polygon
p2.translate(x2,y2);
p3.translate(x2,y2);
drawLine(canvas,x2,y2,p2.x,p2.y,color);
drawLine(canvas,x2,y2,p3.x,p3.y,color);
}
}
/*------------------------------------------------------------------*/
/* Draw a line between two points.
Bresenham's algorithm */
static public void drawLine(double [][]canvas, int x1, int y1,
int x2, int y2, double color) {
int temp;
int dy_neg = 1;
int dx_neg = 1;
int switch_x_y = 0;
int neg_slope = 0;
int tempx, tempy;
int dx = x2 - x1;
if(dx == 0)
if(y1 > y2) {
for(int n = y2; n <= y1; n++) Point(canvas,n,x1,color);
return;
} else {
for(int n = y1; n <= y2; n++) Point(canvas,n,x1,color);
return;
}
int dy = y2 - y1;
if(dy == 0)
if(x1 > x2) {
for(int n = x2; n <= x1; n++) Point(canvas,y1,n,color);
return;
} else {
for(int n = x1; n <= x2; n++) Point(canvas,y1,n,color);
return;
}
float m = (float) dy/dx;
if(m > 1 || m < -1) {
temp = x1;
x1 = y1;
y1 = temp;
temp = x2;
x2 = y2;
y2 = temp;
dx = x2 - x1;
dy = y2 - y1;
m = (float) dy/dx;
switch_x_y = 1;
}
if(x1 > x2) {
temp = x1;
x1 = x2;
x2 = temp;
temp = y1;
y1 = y2;
y2 = temp;
dx = x2 - x1;
dy = y2 - y1;
m = (float) dy/dx;
}
if(m < 0) {
if(dy < 0) {
dy_neg = -1;
dx_neg = 1;
} else {
dy_neg = 1;
dx_neg = -1;
}
neg_slope = 1;
}
int d = 2 * (dy * dy_neg) - (dx * dx_neg);
int incrH = 2 * dy * dy_neg;
int incrHV = 2 * ( (dy * dy_neg) - (dx * dx_neg) );
int x = x1;
int y = y1;
tempx = x;
tempy = y;
if(switch_x_y == 1) {
temp = x;
x = y;
y = temp;
}
Point(canvas,y,x,color);
x = tempx;
y = tempy;
while(x < x2) {
if(d <= 0) {
x++;
d += incrH;
} else {
d += incrHV;
x++;
if(neg_slope == 0) y++;
else y--;
}
tempx = x;
tempy = y;
if (switch_x_y == 1) {
temp = x;
x = y;
y = temp;
}
Point(canvas,y,x,color);
x = tempx;
y = tempy;
}
}
/*------------------------------------------------------------------*/
static public void extractImage(final ImageProcessor ip, double image[]) {
int k=0;
int height=ip.getHeight();
int width =ip.getWidth ();
if (ip instanceof ByteProcessor) {
final byte[] pixels = (byte[])ip.getPixels();
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
image[k] = (double)(pixels[k] & 0xFF);
} else if (ip instanceof ShortProcessor) {
final short[] pixels = (short[])ip.getPixels();
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
if (pixels[k] < (short)0) image[k] = (double)pixels[k] + 65536.0F;
else image[k] = (double)pixels[k];
} else if (ip instanceof FloatProcessor) {
final float[] pixels = (float[])ip.getPixels();
for (int p = 0; p<height*width; p++)
image[p]=pixels[p];
}
}
static public void extractImage(final ImageProcessor ip, double image[][]) {
int k=0;
int height=ip.getHeight();
int width =ip.getWidth ();
if (ip instanceof ByteProcessor) {
final byte[] pixels = (byte[])ip.getPixels();
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
image[y][x] = (double)(pixels[k] & 0xFF);
} else if (ip instanceof ShortProcessor) {
final short[] pixels = (short[])ip.getPixels();
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
if (pixels[k] < (short)0) image[y][x] = (double)pixels[k] + 65536.0F;
else image[y][x] = (double)pixels[k];
} else if (ip instanceof FloatProcessor) {
final float[] pixels = (float[])ip.getPixels();
for (int y = 0; (y < height); y++)
for (int x = 0; (x < width); x++, k++)
image[y][x]=pixels[k];
}
}
/*------------------------------------------------------------------*/
/* Load landmarks from file. */
static public void loadPoints(String filename,
Stack sourceStack, Stack targetStack) {
Point sourcePoint;
Point targetPoint;
try {
final FileReader fr = new FileReader(filename);
final BufferedReader br = new BufferedReader(fr);
String line;
String index;
String xSource;
String ySource;
String xTarget;
String yTarget;
int separatorIndex;
int k = 1;
if (!(line = br.readLine()).equals("Index\txSource\tySource\txTarget\tyTarget")) {
fr.close();
IJ.write("Line " + k + ": 'Index\txSource\tySource\txTarget\tyTarget'");
return;
}
++k;
while ((line = br.readLine()) != null) {
line = line.trim();
separatorIndex = line.indexOf('\t');
if (separatorIndex == -1) {
fr.close();
IJ.write("Line " + k
+ ": #Index# <tab> #xSource# <tab> #ySource# <tab> #xTarget# <tab> #yTarget#");
return;
}
index = line.substring(0, separatorIndex);
index = index.trim();
line = line.substring(separatorIndex);
line = line.trim();
separatorIndex = line.indexOf('\t');
if (separatorIndex == -1) {
fr.close();
IJ.write("Line " + k
+ ": #Index# <tab> #xSource# <tab> #ySource# <tab> #xTarget# <tab> #yTarget#");
return;
}
xSource = line.substring(0, separatorIndex);
xSource = xSource.trim();
line = line.substring(separatorIndex);
line = line.trim();
separatorIndex = line.indexOf('\t');
if (separatorIndex == -1) {
fr.close();
IJ.write("Line " + k
+ ": #Index# <tab> #xSource# <tab> #ySource# <tab> #xTarget# <tab> #yTarget#");
return;
}
ySource = line.substring(0, separatorIndex);
ySource = ySource.trim();
line = line.substring(separatorIndex);
line = line.trim();
separatorIndex = line.indexOf('\t');
if (separatorIndex == -1) {
fr.close();
IJ.write("Line " + k
+ ": #Index# <tab> #xSource# <tab> #ySource# <tab> #xTarget# <tab> #yTarget#");
return;
}
xTarget = line.substring(0, separatorIndex);
xTarget = xTarget.trim();
yTarget = line.substring(separatorIndex);
yTarget = yTarget.trim();
sourcePoint = new Point(Integer.valueOf(xSource).intValue(),
Integer.valueOf(ySource).intValue());
sourceStack.push(sourcePoint);
targetPoint = new Point(Integer.valueOf(xTarget).intValue(),
Integer.valueOf(yTarget).intValue());
targetStack.push(targetPoint);
}
fr.close();
} catch (FileNotFoundException e) {
IJ.error("File not found exception" + e);
return;
} catch (IOException e) {
IJ.error("IOException exception" + e);
return;
} catch (NumberFormatException e) {
IJ.error("Number format exception" + e);
return;
}
}
static public void loadTransformation(String filename,
final double [][]cx, final double [][]cy) {
try {
final FileReader fr = new FileReader(filename);
final BufferedReader br = new BufferedReader(fr);
String line;
// Read number of intervals
line = br.readLine();
int lineN=1;
StringTokenizer st=new StringTokenizer(line,"=");
if (st.countTokens()!=2) {
fr.close();
IJ.write("Line "+lineN+"+: Cannot read number of intervals");
return;
}
st.nextToken();
int intervals=Integer.valueOf(st.nextToken()).intValue();
// Skip next 2 lines
line = br.readLine();
line = br.readLine();
lineN+=2;
// Read the cx coefficients
for (int i= 0; i<intervals+3; i++) {
line = br.readLine(); lineN++;
st=new StringTokenizer(line);
if (st.countTokens()!=intervals+3) {
fr.close();
IJ.write("Line "+lineN+": Cannot read enough coefficients");
return;
}
for (int j=0; j<intervals+3; j++)
cx[i][j]=Double.valueOf(st.nextToken()).doubleValue();
}
// Skip next 2 lines
line = br.readLine();
line = br.readLine();
lineN+=2;
// Read the cy coefficients
for (int i= 0; i<intervals+3; i++) {
line = br.readLine(); lineN++;
st=new StringTokenizer(line);
if (st.countTokens()!=intervals+3) {
fr.close();
IJ.write("Line "+lineN+": Cannot read enough coefficients");
return;
}
for (int j=0; j<intervals+3; j++)
cy[i][j]=Double.valueOf(st.nextToken()).doubleValue();
}
fr.close();
} catch (FileNotFoundException e) {
IJ.error("File not found exception" + e);
return;
} catch (IOException e) {
IJ.error("IOException exception" + e);
return;
} catch (NumberFormatException e) {
IJ.error("Number format exception" + e);
return;
}
}
/*------------------------------------------------------------------*/
static public int numberOfIntervalsOfTransformation(String filename) {
try {
final FileReader fr = new FileReader(filename);
final BufferedReader br = new BufferedReader(fr);
String line;
// Read number of intervals
line = br.readLine();
int lineN=1;
StringTokenizer st=new StringTokenizer(line,"=");
if (st.countTokens()!=2) {
fr.close();
IJ.write("Line "+lineN+"+: Cannot read number of intervals");
return -1;
}
st.nextToken();
int intervals=Integer.valueOf(st.nextToken()).intValue();
fr.close();
return intervals;
} catch (FileNotFoundException e) {
IJ.error("File not found exception" + e);
return -1;
} catch (IOException e) {
IJ.error("IOException exception" + e);
return -1;
} catch (NumberFormatException e) {
IJ.error("Number format exception" + e);
return -1;
}
}
/*------------------------------------------------------------------*/
/* Plot a point in a canvas. */
static public void Point(double [][]canvas, int y, int x, double color) {
if (y<0 || y>=canvas.length) return;
if (x<0 || x>=canvas[0].length) return;
canvas[y][x]=color;
}
/*------------------------------------------------------------------*/
public static void printMatrix(
final String title,
final double [][]array
) {
int Ydim=array.length;
int Xdim=array[0].length;
System.out.println(title);
for (int i=0; i<Ydim; i++) {
for (int j=0; j<Xdim; j++)
System.out.print(array[i][j]+" ");
System.out.println();
}
}
/*------------------------------------------------------------------*/
public static void showImage(
final String title,
final double [] array,
final int Ydim,
final int Xdim
) {
final FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
int ij=0;
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++, ij++)
fp.putPixelValue(j,i,array[ij]);
fp.resetMinAndMax();
final ImagePlus ip=new ImagePlus(title, fp);
final ImageWindow iw=new ImageWindow(ip);
ip.updateImage();
}
/*------------------------------------------------------------------*/
public static void showImage(
final String title,
final double [][]array
) {
int Ydim=array.length;
int Xdim=array[0].length;
final FloatProcessor fp=new FloatProcessor(Xdim,Ydim);
for (int i=0; i<Ydim; i++)
for (int j=0; j<Xdim; j++)
fp.putPixelValue(j,i,array[i][j]);
fp.resetMinAndMax();
final ImagePlus ip=new ImagePlus(title, fp);
final ImageWindow iw=new ImageWindow(ip);
ip.updateImage();
}
} /* End of MiscTools class */
/*====================================================================
| unwarpJPointAction
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJPointAction
extends
ImageCanvas
implements
KeyListener,
MouseListener,
MouseMotionListener
{ /* begin class unwarpJPointAction */
/*....................................................................
Public variables
....................................................................*/
public static final int ADD_CROSS = 0;
public static final int MOVE_CROSS = 1;
public static final int REMOVE_CROSS = 2;
public static final int MASK = 3;
public static final int INVERTMASK = 4;
public static final int FILE = 5;
public static final int STOP = 7;
public static final int MAGNIFIER = 11;
/*....................................................................
Private variables
....................................................................*/
private ImagePlus mainImp;
private ImagePlus secondaryImp;
private unwarpJPointHandler mainPh;
private unwarpJPointHandler secondaryPh;
private unwarpJPointToolbar tb;
private unwarpJDialog dialog;
private long mouseDownTime;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void keyPressed (
final KeyEvent e
) {
if (tb.getCurrentTool()==MASK || tb.getCurrentTool()==INVERTMASK) return;
final Point p = mainPh.getPoint();
if (p == null) return;
final int x = p.x;
final int y = p.y;
switch (e.getKeyCode()) {
case KeyEvent.VK_DELETE:
case KeyEvent.VK_BACK_SPACE:
mainPh.removePoint();
secondaryPh.removePoint();
updateAndDraw();
break;
case KeyEvent.VK_DOWN:
mainPh.movePoint(mainImp.getWindow().getCanvas().screenX(x),
mainImp.getWindow().getCanvas().screenY(y
+ (int)Math.ceil(1.0 / mainImp.getWindow().getCanvas().getMagnification())));
mainImp.setRoi(mainPh);
break;
case KeyEvent.VK_LEFT:
mainPh.movePoint(mainImp.getWindow().getCanvas().screenX(x
- (int)Math.ceil(1.0 / mainImp.getWindow().getCanvas().getMagnification())),
mainImp.getWindow().getCanvas().screenY(y));
mainImp.setRoi(mainPh);
break;
case KeyEvent.VK_RIGHT:
mainPh.movePoint(mainImp.getWindow().getCanvas().screenX(x
+ (int)Math.ceil(1.0 / mainImp.getWindow().getCanvas().getMagnification())),
mainImp.getWindow().getCanvas().screenY(y));
mainImp.setRoi(mainPh);
break;
case KeyEvent.VK_TAB:
mainPh.nextPoint();
secondaryPh.nextPoint();
updateAndDraw();
break;
case KeyEvent.VK_UP:
mainPh.movePoint(mainImp.getWindow().getCanvas().screenX(x),
mainImp.getWindow().getCanvas().screenY(y
- (int)Math.ceil(1.0 / mainImp.getWindow().getCanvas().getMagnification())));
mainImp.setRoi(mainPh);
break;
}
} /* end keyPressed */
/*------------------------------------------------------------------*/
public void keyReleased (
final KeyEvent e
) {
} /* end keyReleased */
/*------------------------------------------------------------------*/
public void keyTyped (
final KeyEvent e
) {
} /* end keyTyped */
/*------------------------------------------------------------------*/
public void mouseClicked (
final MouseEvent e
) {
} /* end mouseClicked */
/*------------------------------------------------------------------*/
public void mouseDragged (
final MouseEvent e
) {
final int x = e.getX();
final int y = e.getY();
if (tb.getCurrentTool() == MOVE_CROSS) {
mainPh.movePoint(x, y);
updateAndDraw();
}
mouseMoved(e);
} /* end mouseDragged */
/*------------------------------------------------------------------*/
public void mouseEntered (
final MouseEvent e
) {
WindowManager.setCurrentWindow(mainImp.getWindow());
mainImp.getWindow().toFront();
updateAndDraw();
} /* end mouseEntered */
/*------------------------------------------------------------------*/
public void mouseExited (
final MouseEvent e
) {
IJ.showStatus("");
} /* end mouseExited */
/*------------------------------------------------------------------*/
public void mouseMoved (
final MouseEvent e
) {
setControl();
final int x = mainImp.getWindow().getCanvas().offScreenX(e.getX());
final int y = mainImp.getWindow().getCanvas().offScreenY(e.getY());
IJ.showStatus(mainImp.getLocationAsString(x, y) + getValueAsString(x, y));
} /* end mouseMoved */
/*------------------------------------------------------------------*/
public void mousePressed (
final MouseEvent e
) {
if (dialog.isFinalActionLaunched()) return;
int x = e.getX(),xp;
int y = e.getY(),yp;
int currentPoint;
boolean doubleClick = (System.currentTimeMillis() - mouseDownTime) <= 250L;
mouseDownTime = System.currentTimeMillis();
switch (tb.getCurrentTool()) {
case ADD_CROSS:
xp=mainImp.getWindow().getCanvas().offScreenX(x);
yp=mainImp.getWindow().getCanvas().offScreenY(y);
mainPh.addPoint(xp,yp);
xp = positionX(mainImp, secondaryImp, mainImp.getWindow().getCanvas().offScreenX(x));
yp = positionY(mainImp, secondaryImp, mainImp.getWindow().getCanvas().offScreenY(y));
secondaryPh.addPoint(xp, yp);
updateAndDraw();
break;
case MOVE_CROSS:
currentPoint = mainPh.findClosest(x, y);
secondaryPh.setCurrentPoint(currentPoint);
updateAndDraw();
break;
case REMOVE_CROSS:
currentPoint = mainPh.findClosest(x, y);
mainPh.removePoint(currentPoint);
secondaryPh.removePoint(currentPoint);
updateAndDraw();
break;
case MASK:
case INVERTMASK:
if (mainPh.canAddMaskPoints()) {
if (!doubleClick) {
if (dialog.isClearMaskSet()) {
mainPh.clearMask();
dialog.setClearMask(false);
dialog.ungrayImage(this);
}
x = positionX(mainImp, secondaryImp, mainImp.getWindow().getCanvas().offScreenX(x));
y = positionY(mainImp, secondaryImp, mainImp.getWindow().getCanvas().offScreenY(y));
mainPh.addMaskPoint(x,y);
} else
mainPh.closeMask(tb.getCurrentTool());
updateAndDraw();
} else {
IJ.error("A mask cannot be manually assigned since the mask was already in the stack");
}
break;
case MAGNIFIER:
final int flags = e.getModifiers();
if ((flags & (Event.ALT_MASK | Event.META_MASK | Event.CTRL_MASK)) != 0) {
mainImp.getWindow().getCanvas().zoomOut(x, y);
}
else {
mainImp.getWindow().getCanvas().zoomIn(x, y);
}
break;
}
} /* end mousePressed */
/*------------------------------------------------------------------*/
public void mouseReleased (
final MouseEvent e
) {
} /* end mouseReleased */
/*------------------------------------------------------------------*/
public void setSecondaryPointHandler (
final ImagePlus secondaryImp,
final unwarpJPointHandler secondaryPh
) {
this.secondaryImp = secondaryImp;
this.secondaryPh = secondaryPh;
} /* end setSecondaryPointHandler */
/*------------------------------------------------------------------*/
public unwarpJPointAction (
final ImagePlus imp,
final unwarpJPointHandler ph,
final unwarpJPointToolbar tb,
final unwarpJDialog dialog
) {
super(imp);
this.mainImp = imp;
this.mainPh = ph;
this.tb = tb;
this.dialog = dialog;
} /* end unwarpJPointAction */
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private String getValueAsString (
final int x,
final int y
) {
final Calibration cal = mainImp.getCalibration();
final int[] v = mainImp.getPixel(x, y);
final int mainImptype=mainImp.getType();
if (mainImptype==mainImp.GRAY8 || mainImptype==mainImp.GRAY16) {
final double cValue = cal.getCValue(v[0]);
if (cValue==v[0]) {
return(", value=" + v[0]);
}
else {
return(", value=" + IJ.d2s(cValue) + " (" + v[0] + ")");
}
} else if (mainImptype==mainImp.GRAY32) {
return(", value=" + Float.intBitsToFloat(v[0]));
} else if (mainImptype==mainImp.COLOR_256) {
return(", index=" + v[3] + ", value=" + v[0] + "," + v[1] + "," + v[2]);
} else if (mainImptype==mainImp.COLOR_RGB) {
return(", value=" + v[0] + "," + v[1] + "," + v[2]);
} else {
return("");
}
} /* end getValueAsString */
/*------------------------------------------------------------------*/
private int positionX (
final ImagePlus imp1,
final ImagePlus imp2,
final int x
) {
return((x * imp2.getWidth()) / imp1.getWidth());
} /* end PositionX */
/*------------------------------------------------------------------*/
private int positionY (
final ImagePlus imp1,
final ImagePlus imp2,
final int y
) {
return((y * imp2.getHeight()) / imp1.getHeight());
} /* end PositionY */
/*------------------------------------------------------------------*/
private void setControl (
) {
switch (tb.getCurrentTool()) {
case ADD_CROSS:
mainImp.getWindow().getCanvas().setCursor(crosshairCursor);
break;
case FILE:
case MAGNIFIER:
case MOVE_CROSS:
case REMOVE_CROSS:
case MASK:
case INVERTMASK:
case STOP:
mainImp.getWindow().getCanvas().setCursor(defaultCursor);
break;
}
} /* end setControl */
/*------------------------------------------------------------------*/
private void updateAndDraw (
) {
mainImp.setRoi(mainPh);
secondaryImp.setRoi(secondaryPh);
} /* end updateAndDraw */
} /* end class unwarpJPointAction */
/*====================================================================
| unwarpJPointHandler
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJPointHandler
extends
Roi
{ /* begin class unwarpJPointHandler */
/*....................................................................
Private variables
....................................................................*/
private static final int CROSS_HALFSIZE = 5;
// Colors
private static final int GAMUT = 1024;
private final Color spectrum[] = new Color[GAMUT];
private final boolean usedColor[] = new boolean[GAMUT];
private final Vector listColors = new Vector(0, 16);
private int currentColor = 0;
// List of crosses
private final Vector listPoints = new Vector(0, 16);
private int currentPoint = -1;
private int numPoints = 0;
private boolean started = false;
// List of points for the mask
private final Vector listMaskPoints = new Vector(0,16);
private boolean maskClosed = false;
// Some useful references
private ImagePlus imp;
private unwarpJPointAction pa;
private unwarpJPointToolbar tb;
private unwarpJMask mask;
private unwarpJDialog dialog;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void addMaskPoint (
final int x,
final int y
) {
if (maskClosed) return;
final Point p = new Point(x, y);
listMaskPoints.addElement(p);
}
/*------------------------------------------------------------------*/
public void addPoint (
final int x,
final int y
) {
if (numPoints < GAMUT) {
final Point p = new Point(x, y);
listPoints.addElement(p);
if (!usedColor[currentColor]) {
usedColor[currentColor] = true;
}
else {
int k;
for (k = 0; (k < GAMUT); k++) {
currentColor++;
currentColor &= GAMUT - 1;
if (!usedColor[currentColor]) {
break;
}
}
if (GAMUT <= k) {
throw new IllegalStateException("Unexpected lack of available colors");
}
}
int stirredColor = 0;
int c = currentColor;
for (int k = 0; (k < (int)Math.round(Math.log((double)GAMUT) / Math.log(2.0))); k++) {
stirredColor <<= 1;
stirredColor |= (c & 1);
c >>= 1;
}
listColors.addElement(new Integer(stirredColor));
currentColor++;
currentColor &= GAMUT - 1;
currentPoint = numPoints;
numPoints++;
}
else {
IJ.error("Maximum number of points reached");
}
} /* end addPoint */
/*------------------------------------------------------------------*/
/* False if the image is coming from a stack */
public boolean canAddMaskPoints() {
return !mask.isFromStack();
}
/*------------------------------------------------------------------*/
public void clearMask () {
// Clear mask information in this object
listMaskPoints.removeAllElements();
maskClosed=false;
mask.clearMask();
}
/*------------------------------------------------------------------*/
public void closeMask (int tool) {
listMaskPoints.addElement(listMaskPoints.elementAt(0));
maskClosed=true;
mask.setMaskPoints(listMaskPoints);
mask.fillMask(tool);
dialog.grayImage(this);
}
/*------------------------------------------------------------------*/
public void draw (
final Graphics g
) {
// Draw landmarks
if (started) {
final double mag = (double)ic.getMagnification();
final int dx = (int)(mag / 2.0);
final int dy = (int)(mag / 2.0);
for (int k = 0; (k < numPoints); k++) {
final Point p = (Point)listPoints.elementAt(k);
g.setColor(spectrum[((Integer)listColors.elementAt(k)).intValue()]);
if (k == currentPoint) {
if (WindowManager.getCurrentImage() == imp) {
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y - 1) + dy,
ic.screenX(p.x - 1) + dx,
ic.screenY(p.y - 1) + dy);
g.drawLine(ic.screenX(p.x - 1) + dx,
ic.screenY(p.y - 1) + dy,
ic.screenX(p.x - 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE - 1) + dy);
g.drawLine(ic.screenX(p.x - 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE - 1) + dy,
ic.screenX(p.x + 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE - 1) + dy);
g.drawLine(ic.screenX(p.x + 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE - 1) + dy,
ic.screenX(p.x + 1) + dx,
ic.screenY(p.y - 1) + dy);
g.drawLine(ic.screenX(p.x + 1) + dx,
ic.screenY(p.y - 1) + dy,
ic.screenX(p.x + CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y - 1) + dy);
g.drawLine(ic.screenX(p.x + CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y - 1) + dy,
ic.screenX(p.x + CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y + 1) + dy);
g.drawLine(ic.screenX(p.x + CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y + 1) + dy,
ic.screenX(p.x + 1) + dx,
ic.screenY(p.y + 1) + dy);
g.drawLine(ic.screenX(p.x + 1) + dx,
ic.screenY(p.y + 1) + dy,
ic.screenX(p.x + 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE + 1) + dy);
g.drawLine(ic.screenX(p.x + 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE + 1) + dy,
ic.screenX(p.x - 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE + 1) + dy);
g.drawLine(ic.screenX(p.x - 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE + 1) + dy,
ic.screenX(p.x - 1) + dx,
ic.screenY(p.y + 1) + dy);
g.drawLine(ic.screenX(p.x - 1) + dx,
ic.screenY(p.y + 1) + dy,
ic.screenX(p.x - CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y + 1) + dy);
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y + 1) + dy,
ic.screenX(p.x - CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y - 1) + dy);
if (1.0 < ic.getMagnification()) {
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE) + dx,
ic.screenY(p.y) + dy,
ic.screenX(p.x + CROSS_HALFSIZE) + dx,
ic.screenY(p.y) + dy);
g.drawLine(ic.screenX(p.x) + dx,
ic.screenY(p.y - CROSS_HALFSIZE) + dy,
ic.screenX(p.x) + dx,
ic.screenY(p.y + CROSS_HALFSIZE) + dy);
}
}
else {
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE + 1) + dy,
ic.screenX(p.x + CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE - 1) + dy);
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE + 1) + dx,
ic.screenY(p.y + CROSS_HALFSIZE - 1) + dy,
ic.screenX(p.x + CROSS_HALFSIZE - 1) + dx,
ic.screenY(p.y - CROSS_HALFSIZE + 1) + dy);
}
}
else {
g.drawLine(ic.screenX(p.x - CROSS_HALFSIZE) + dx,
ic.screenY(p.y) + dy,
ic.screenX(p.x + CROSS_HALFSIZE) + dx,
ic.screenY(p.y) + dy);
g.drawLine(ic.screenX(p.x) + dx,
ic.screenY(p.y - CROSS_HALFSIZE) + dy,
ic.screenX(p.x) + dx,
ic.screenY(p.y + CROSS_HALFSIZE) + dy);
}
}
if (updateFullWindow) {
updateFullWindow = false;
imp.draw();
}
}
// Draw mask
int numberMaskPoints=listMaskPoints.size();
if (numberMaskPoints!=0) {
final double mag = (double)ic.getMagnification();
final int dx = (int)(mag / 2.0);
final int dy = (int)(mag / 2.0);
int CIRCLE_RADIUS=CROSS_HALFSIZE/2;
int CIRCLE_DIAMETER=2*CIRCLE_RADIUS;
for (int i=0; i<numberMaskPoints; i++) {
final Point p = (Point)listMaskPoints.elementAt(i);
g.setColor(Color.yellow);
g.drawOval(ic.screenX(p.x)-CIRCLE_RADIUS+dx, ic.screenY(p.y)-CIRCLE_RADIUS+dy,
CIRCLE_DIAMETER, CIRCLE_DIAMETER);
if (i!=0) {
Point previous_p=(Point)listMaskPoints.elementAt(i-1);
g.drawLine(ic.screenX(p.x)+dx,ic.screenY(p.y)+dy,
ic.screenX(previous_p.x)+dx,ic.screenY(previous_p.y)+dy);
}
}
}
} /* end draw */
/*------------------------------------------------------------------*/
public int findClosest (
int x,
int y
) {
if (numPoints == 0) {
return(currentPoint);
}
x = ic.offScreenX(x);
y = ic.offScreenY(y);
Point p = new Point((Point)listPoints.elementAt(currentPoint));
double distance = (double)(x - p.x) * (double)(x - p.x)
+ (double)(y - p.y) * (double)(y - p.y);
for (int k = 0; (k < numPoints); k++) {
p = (Point)listPoints.elementAt(k);
final double candidate = (double)(x - p.x) * (double)(x - p.x)
+ (double)(y - p.y) * (double)(y - p.y);
if (candidate < distance) {
distance = candidate;
currentPoint = k;
}
}
return(currentPoint);
} /* end findClosest */
/*------------------------------------------------------------------*/
public Point getPoint (
) {
return((0 <= currentPoint) ? (Point)listPoints.elementAt(currentPoint) : (null));
} /* end getPoint */
/*------------------------------------------------------------------*/
public unwarpJPointAction getPointAction () {return pa;}
/*------------------------------------------------------------------*/
public int getCurrentPoint (
) {
return(currentPoint);
} /* end getCurrentPoint */
/*------------------------------------------------------------------*/
public Vector getPoints (
) {
return(listPoints);
} /* end getPoints */
/*------------------------------------------------------------------*/
public void killListeners (
) {
final ImageWindow iw = imp.getWindow();
final ImageCanvas ic = iw.getCanvas();
ic.removeKeyListener(pa);
ic.removeMouseListener(pa);
ic.removeMouseMotionListener(pa);
ic.addMouseMotionListener(ic);
ic.addMouseListener(ic);
ic.addKeyListener(IJ.getInstance());
} /* end killListeners */
/*------------------------------------------------------------------*/
public void movePoint (
int x,
int y
) {
if (0 <= currentPoint) {
x = ic.offScreenX(x);
y = ic.offScreenY(y);
x = (x < 0) ? (0) : (x);
x = (imp.getWidth() <= x) ? (imp.getWidth() - 1) : (x);
y = (y < 0) ? (0) : (y);
y = (imp.getHeight() <= y) ? (imp.getHeight() - 1) : (y);
listPoints.removeElementAt(currentPoint);
final Point p = new Point(x, y);
listPoints.insertElementAt(p, currentPoint);
}
} /* end movePoint */
/*------------------------------------------------------------------*/
public void nextPoint (
) {
currentPoint = (currentPoint == (numPoints - 1)) ? (0) : (currentPoint + 1);
} /* end nextPoint */
/*------------------------------------------------------------------*/
public void removePoint (
) {
if (0 < numPoints) {
listPoints.removeElementAt(currentPoint);
usedColor[((Integer)listColors.elementAt(currentPoint)).intValue()] = false;
listColors.removeElementAt(currentPoint);
numPoints--;
}
currentPoint = numPoints - 1;
if (currentPoint < 0) {
tb.setTool(pa.ADD_CROSS);
}
} /* end removePoint */
/*------------------------------------------------------------------*/
public void removePoint (
final int k
) {
if (0 < numPoints) {
listPoints.removeElementAt(k);
usedColor[((Integer)listColors.elementAt(k)).intValue()] = false;
listColors.removeElementAt(k);
numPoints--;
}
currentPoint = numPoints - 1;
if (currentPoint < 0) {
tb.setTool(pa.ADD_CROSS);
}
} /* end removePoint */
/*------------------------------------------------------------------*/
public void removePoints (
) {
listPoints.removeAllElements();
listColors.removeAllElements();
for (int k = 0; (k < GAMUT); k++) {
usedColor[k] = false;
}
currentColor = 0;
numPoints = 0;
currentPoint = -1;
tb.setTool(pa.ADD_CROSS);
imp.setRoi(this);
} /* end removePoints */
/*------------------------------------------------------------------*/
public void setCurrentPoint (
final int currentPoint
) {
this.currentPoint = currentPoint;
} /* end setCurrentPoint */
/*------------------------------------------------------------------*/
public void setTestSourceSet (
final int set
) {
removePoints();
switch(set) {
case 1: // Deformed_Lena 1
addPoint(11,11);
addPoint(200,6);
addPoint(197,204);
addPoint(121,111);
break;
case 2: // Deformed_Lena 1
addPoint(6,6);
addPoint(202,7);
addPoint(196,210);
addPoint(10,214);
addPoint(120,112);
addPoint(68,20);
addPoint(63,163);
addPoint(186,68);
break;
}
} /* end setTestset */
/*------------------------------------------------------------------*/
public void setTestTargetSet (
final int set
) {
removePoints();
switch(set) {
case 1:
addPoint(11,11);
addPoint(185,15);
addPoint(154,200);
addPoint(123,92);
break;
case 2: // Deformed_Lena 1
addPoint(6,6);
addPoint(185,14);
addPoint(154,200);
addPoint(3,178);
addPoint(121,93);
addPoint(67,14);
addPoint(52,141);
addPoint(178,68);
break;
}
} /* end setTestset */
/*------------------------------------------------------------------*/
public void setSecondaryPointHandler (
final ImagePlus secondaryImp,
final unwarpJPointHandler secondaryPh
) {
pa.setSecondaryPointHandler(secondaryImp, secondaryPh);
} /* end setSecondaryPointHandler */
/*------------------------------------------------------------------*/
/* Constructor with graphical capabilities */
public unwarpJPointHandler (
final ImagePlus imp,
final unwarpJPointToolbar tb,
final unwarpJMask mask,
final unwarpJDialog dialog
) {
super(0, 0, imp.getWidth(), imp.getHeight(), imp);
this.imp = imp;
this.tb = tb;
this.dialog=dialog;
pa = new unwarpJPointAction(imp, this, tb, dialog);
final ImageWindow iw = imp.getWindow();
final ImageCanvas ic = iw.getCanvas();
iw.requestFocus();
iw.removeKeyListener(IJ.getInstance());
iw.addKeyListener(pa);
ic.removeMouseMotionListener(ic);
ic.removeMouseListener(ic);
ic.removeKeyListener(IJ.getInstance());
ic.addKeyListener(pa);
ic.addMouseListener(pa);
ic.addMouseMotionListener(pa);
setSpectrum();
started = true;
this.mask=mask;
clearMask();
} /* end unwarpJPointHandler */
/* Constructor without graphical capabilities */
public unwarpJPointHandler (
final ImagePlus imp
) {
super(0, 0, imp.getWidth(), imp.getHeight(), imp);
this.imp = imp;
tb = null;
dialog=null;
pa = null;
started = true;
mask=null;
} /* end unwarpJPointHandler */
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void setSpectrum (
) {
final int bound1 = GAMUT / 6;
final int bound2 = GAMUT / 3;
final int bound3 = GAMUT / 2;
final int bound4 = (2 * GAMUT) / 3;
final int bound5 = (5 * GAMUT) / 6;
final int bound6 = GAMUT;
final float gamutChunk1 = (float)bound1;
final float gamutChunk2 = (float)(bound2 - bound1);
final float gamutChunk3 = (float)(bound3 - bound2);
final float gamutChunk4 = (float)(bound4 - bound3);
final float gamutChunk5 = (float)(bound5 - bound4);
final float gamutChunk6 = (float)(bound6 - bound5);
int k = 0;
do {
spectrum[k] = new Color(1.0F, (float)k / gamutChunk1, 0.0F);
usedColor[k] = false;
} while (++k < bound1);
do {
spectrum[k] = new Color(1.0F - (float)(k - bound1) / gamutChunk2, 1.0F, 0.0F);
usedColor[k] = false;
} while (++k < bound2);
do {
spectrum[k] = new Color(0.0F, 1.0F, (float)(k - bound2) / gamutChunk3);
usedColor[k] = false;
} while (++k < bound3);
do {
spectrum[k] = new Color(0.0F, 1.0F - (float)(k - bound3) / gamutChunk4, 1.0F);
usedColor[k] = false;
} while (++k < bound4);
do {
spectrum[k] = new Color((float)(k - bound4) / gamutChunk5, 0.0F, 1.0F);
usedColor[k] = false;
} while (++k < bound5);
do {
spectrum[k] = new Color(1.0F, 0.0F, 1.0F - (float)(k - bound5) / gamutChunk6);
usedColor[k] = false;
} while (++k < bound6);
} /* end setSpectrum */
} /* end class unwarpJPointHandler */
/*====================================================================
| unwarpJPointToolbar
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJPointToolbar
extends
Canvas
implements
MouseListener
{ /* begin class unwarpJPointToolbar */
/*....................................................................
Private variables
....................................................................*/
private static final int NUM_TOOLS = 19;
private static final int SIZE = 22;
private static final int OFFSET = 3;
private static final Color gray = Color.lightGray;
private static final Color brighter = gray.brighter();
private static final Color darker = gray.darker();
private static final Color evenDarker = darker.darker();
private final boolean[] down = new boolean[NUM_TOOLS];
private Graphics g;
private ImagePlus sourceImp;
private ImagePlus targetImp;
private Toolbar previousInstance;
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
private unwarpJPointToolbar instance;
private long mouseDownTime;
private int currentTool = unwarpJPointAction.ADD_CROSS;
private int x;
private int y;
private int xOffset;
private int yOffset;
private unwarpJDialog dialog;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public int getCurrentTool (
) {
return(currentTool);
} /* getCurrentTool */
/*------------------------------------------------------------------*/
public void mouseClicked (
final MouseEvent e
) {
} /* end mouseClicked */
/*------------------------------------------------------------------*/
public void mouseEntered (
final MouseEvent e
) {
} /* end mouseEntered */
/*------------------------------------------------------------------*/
public void mouseExited (
final MouseEvent e
) {
} /* end mouseExited */
/*------------------------------------------------------------------*/
public void mousePressed (
final MouseEvent e
) {
final int x = e.getX();
final int y = e.getY();
int newTool = 0;
for (int i = 0; (i < NUM_TOOLS); i++) {
if (((i * SIZE) < x) && (x < (i * SIZE + SIZE))) {
newTool = i;
}
}
boolean doubleClick = ((newTool == getCurrentTool())
&& ((System.currentTimeMillis() - mouseDownTime) <= 500L)
&& (newTool == unwarpJPointAction.REMOVE_CROSS));
mouseDownTime = System.currentTimeMillis();
if (newTool==unwarpJPointAction.STOP && !dialog.isFinalActionLaunched())
return;
if (newTool!=unwarpJPointAction.STOP && dialog.isFinalActionLaunched())
return;
setTool(newTool);
if (doubleClick) {
unwarpJClearAll clearAllDialog = new unwarpJClearAll(IJ.getInstance(),
sourceImp, targetImp, sourcePh, targetPh);
GUI.center(clearAllDialog);
clearAllDialog.setVisible(true);
setTool(unwarpJPointAction.ADD_CROSS);
clearAllDialog.dispose();
}
switch (newTool) {
case unwarpJPointAction.FILE:
unwarpJFile fileDialog = new unwarpJFile(IJ.getInstance(),
sourceImp, targetImp, sourcePh, targetPh,dialog);
GUI.center(fileDialog);
fileDialog.setVisible(true);
setTool(unwarpJPointAction.ADD_CROSS);
fileDialog.dispose();
break;
case unwarpJPointAction.MASK:
case unwarpJPointAction.INVERTMASK:
dialog.setClearMask(true);
break;
case unwarpJPointAction.STOP:
dialog.setStopRegistration();
break;
}
} /* mousePressed */
/*------------------------------------------------------------------*/
public void mouseReleased (
final MouseEvent e
) {
} /* end mouseReleased */
/*------------------------------------------------------------------*/
public void paint (
final Graphics g
) {
for (int i = 0; (i < NUM_TOOLS); i++) {
drawButton(g, i);
}
} /* paint */
/*------------------------------------------------------------------*/
public void restorePreviousToolbar (
) {
final Container container = instance.getParent();
final Component[] component = container.getComponents();
for (int i = 0; (i < component.length); i++) {
if (component[i] == instance) {
container.remove(instance);
container.add(previousInstance, i);
container.validate();
break;
}
}
} /* end restorePreviousToolbar */
/*------------------------------------------------------------------*/
public void setAllUp () {
for (int i=0; i<NUM_TOOLS; i++) down[i]=false;
}
/*------------------------------------------------------------------*/
public void setSource (
final ImagePlus sourceImp,
final unwarpJPointHandler sourcePh
) {
this.sourceImp = sourceImp;
this.sourcePh = sourcePh;
} /* end setSource */
/*------------------------------------------------------------------*/
public void setTarget (
final ImagePlus targetImp,
final unwarpJPointHandler targetPh
) {
this.targetImp = targetImp;
this.targetPh = targetPh;
} /* end setTarget */
/*------------------------------------------------------------------*/
public void setTool (
final int tool
) {
if (tool == currentTool) {
return;
}
down[tool] = true;
down[currentTool] = false;
Graphics g = this.getGraphics();
drawButton(g, currentTool);
drawButton(g, tool);
g.dispose();
showMessage(tool);
currentTool = tool;
} /* end setTool */
/*------------------------------------------------------------------*/
public unwarpJPointToolbar (
final Toolbar previousToolbar,
final unwarpJDialog dialog
) {
previousInstance = previousToolbar;
this.dialog = dialog;
instance = this;
final Container container = previousToolbar.getParent();
final Component[] component = container.getComponents();
for (int i = 0; (i < component.length); i++) {
if (component[i] == previousToolbar) {
container.remove(previousToolbar);
container.add(this, i);
break;
}
}
resetButtons();
down[currentTool] = true;
setTool(currentTool);
setForeground(evenDarker);
setBackground(gray);
addMouseListener(this);
container.validate();
} /* end unwarpJPointToolbar */
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void d (
int x,
int y
) {
x += xOffset;
y += yOffset;
g.drawLine(this.x, this.y, x, y);
this.x = x;
this.y = y;
} /* end d */
/*------------------------------------------------------------------*/
private void drawButton (
final Graphics g,
final int tool
) {
fill3DRect(g, tool * SIZE + 1, 1, SIZE, SIZE - 1, !down[tool]);
if (tool==unwarpJPointAction.STOP && !dialog.isFinalActionLaunched())
return;
if (tool!=unwarpJPointAction.STOP && dialog.isFinalActionLaunched())
return;
g.setColor(Color.black);
int x = tool * SIZE + OFFSET;
int y = OFFSET;
if (down[tool]) {
x++;
y++;
}
this.g = g;
// Plygon for the mask
int px[]=new int[5]; px[0]=x+4;px[1]=x+ 4;px[2]=x+14;px[3]=x+ 9;px[4]=x+14;
int py[]=new int[5]; py[0]=y+3;py[1]=y+13;py[2]=y+13;py[3]=y+ 8;py[4]=y+ 3;
switch (tool) {
case unwarpJPointAction.ADD_CROSS:
xOffset = x;
yOffset = y;
m(7, 0);
d(7, 1);
m(6, 2);
d(6, 3);
m(8, 2);
d(8, 3);
m(5, 4);
d(5, 5);
m(9, 4);
d(9, 5);
m(4, 6);
d(4, 8);
m(10, 6);
d(10, 8);
m(5, 9);
d(5, 14);
m(9, 9);
d(9, 14);
m(7, 4);
d(7, 6);
m(7, 8);
d(7, 8);
m(4, 11);
d(10, 11);
g.fillRect(x + 6, y + 12, 3, 3);
m(11, 13);
d(15, 13);
m(13, 11);
d(13, 15);
break;
case unwarpJPointAction.FILE:
xOffset = x;
yOffset = y;
m(3, 1);
d(9, 1);
d(9, 4);
d(12, 4);
d(12, 14);
d(3, 14);
d(3, 1);
m(10, 2);
d(11, 3);
m(5, 4);
d(7, 4);
m(5, 6);
d(10, 6);
m(5, 8);
d(10, 8);
m(5, 10);
d(10, 10);
m(5, 12);
d(10, 12);
break;
case unwarpJPointAction.MAGNIFIER:
xOffset = x + 2;
yOffset = y + 2;
m(3, 0);
d(3, 0);
d(5, 0);
d(8, 3);
d(8, 5);
d(7, 6);
d(7, 7);
d(6, 7);
d(5, 8);
d(3, 8);
d(0, 5);
d(0, 3);
d(3, 0);
m(8, 8);
d(9, 8);
d(13, 12);
d(13, 13);
d(12, 13);
d(8, 9);
d(8, 8);
break;
case unwarpJPointAction.MOVE_CROSS:
xOffset = x;
yOffset = y;
m(1, 1);
d(1, 10);
m(2, 2);
d(2, 9);
m(3, 3);
d(3, 8);
m(4, 4);
d(4, 7);
m(5, 5);
d(5, 7);
m(6, 6);
d(6, 7);
m(7, 7);
d(7, 7);
m(11, 5);
d(11, 6);
m(10, 7);
d(10, 8);
m(12, 7);
d(12, 8);
m(9, 9);
d(9, 11);
m(13, 9);
d(13, 11);
m(10, 12);
d(10, 15);
m(12, 12);
d(12, 15);
m(11, 9);
d(11, 10);
m(11, 13);
d(11, 15);
m(9, 13);
d(13, 13);
break;
case unwarpJPointAction.REMOVE_CROSS:
xOffset = x;
yOffset = y;
m(7, 0);
d(7, 1);
m(6, 2);
d(6, 3);
m(8, 2);
d(8, 3);
m(5, 4);
d(5, 5);
m(9, 4);
d(9, 5);
m(4, 6);
d(4, 8);
m(10, 6);
d(10, 8);
m(5, 9);
d(5, 14);
m(9, 9);
d(9, 14);
m(7, 4);
d(7, 6);
m(7, 8);
d(7, 8);
m(4, 11);
d(10, 11);
g.fillRect(x + 6, y + 12, 3, 3);
m(11, 13);
d(15, 13);
break;
case unwarpJPointAction.MASK:
xOffset = x;
yOffset = y;
g.fillPolygon(px,py,5);
break;
case unwarpJPointAction.INVERTMASK:
xOffset = x;
yOffset = y;
g.fillRect(x + 1, y + 1, 15, 15);
g.setColor(gray);
g.fillPolygon(px,py,5);
g.setColor(Color.black);
break;
case unwarpJPointAction.STOP:
xOffset = x;
yOffset = y;
// Octogon
m( 1, 5);
d( 1, 11);
d( 5, 15);
d(11, 15);
d(15, 11);
d(15, 5);
d(11, 1);
d( 5, 1);
d( 1, 5);
// S
m( 5, 6);
d( 3, 6);
d( 3, 8);
d( 5, 8);
d( 5, 10);
d( 3, 10);
// T
m( 6, 6);
d( 6, 8);
m( 7, 6);
d( 7, 10);
// O
m(11, 6);
d( 9, 6);
d( 9, 10);
d(11, 10);
d(11, 6);
// P
m(12, 10);
d(12, 6);
d(14, 6);
d(14, 8);
d(12, 8);
break;
}
} /* end drawButton */
/*------------------------------------------------------------------*/
private void fill3DRect (
final Graphics g,
final int x,
final int y,
final int width,
final int height,
final boolean raised
) {
if (raised) {
g.setColor(gray);
}
else {
g.setColor(darker);
}
g.fillRect(x + 1, y + 1, width - 2, height - 2);
g.setColor((raised) ? (brighter) : (evenDarker));
g.drawLine(x, y, x, y + height - 1);
g.drawLine(x + 1, y, x + width - 2, y);
g.setColor((raised) ? (evenDarker) : (brighter));
g.drawLine(x + 1, y + height - 1, x + width - 1, y + height - 1);
g.drawLine(x + width - 1, y, x + width - 1, y + height - 2);
} /* end fill3DRect */
/*------------------------------------------------------------------*/
private void m (
final int x,
final int y
) {
this.x = xOffset + x;
this.y = yOffset + y;
} /* end m */
/*------------------------------------------------------------------*/
private void resetButtons (
) {
for (int i = 0; (i < NUM_TOOLS); i++) {
down[i] = false;
}
} /* end resetButtons */
/*------------------------------------------------------------------*/
private void showMessage (
final int tool
) {
switch (tool) {
case unwarpJPointAction.ADD_CROSS:
IJ.showStatus("Add crosses");
return;
case unwarpJPointAction.FILE:
IJ.showStatus("Export/import list of points");
return;
case unwarpJPointAction.MAGNIFIER:
IJ.showStatus("Magnifying glass");
return;
case unwarpJPointAction.MOVE_CROSS:
IJ.showStatus("Move crosses");
return;
case unwarpJPointAction.REMOVE_CROSS:
IJ.showStatus("Remove crosses");
return;
case unwarpJPointAction.MASK:
IJ.showStatus("Draw a mask");
return;
case unwarpJPointAction.STOP:
IJ.showStatus("Stop registration");
return;
default:
IJ.showStatus("Undefined operation");
return;
}
} /* end showMessage */
} /* end class unwarpJPointToolbar */
/*====================================================================
| unwarpJProgressBar
\===================================================================*/
/*********************************************************************
* This class implements the interactions when dealing with ImageJ's
* progress bar.
********************************************************************/
class unwarpJProgressBar
{ /* begin class unwarpJProgressBar */
/*....................................................................
Private variables
....................................................................*/
/*********************************************************************
* Same time constant than in ImageJ version 1.22
********************************************************************/
private static final long TIME_QUANTUM = 50L;
private static volatile long lastTime = System.currentTimeMillis();
private static volatile int completed = 0;
private static volatile int workload = 0;
/*....................................................................
Public methods
....................................................................*/
/*********************************************************************
* Extend the amount of work to perform by <code>batch</code>.
*
* @param batch Additional amount of work that need be performed.
********************************************************************/
public static synchronized void addWorkload (
final int batch
) {
workload += batch;
} /* end addWorkload */
/*********************************************************************
* Erase the progress bar and cancel pending operations.
********************************************************************/
public static synchronized void resetProgressBar (
) {
final long timeStamp = System.currentTimeMillis();
if ((timeStamp - lastTime) < TIME_QUANTUM) {
try {
Thread.sleep(TIME_QUANTUM - timeStamp + lastTime);
} catch (InterruptedException e) {
IJ.error("Unexpected interruption exception" + e);
}
}
lastTime = timeStamp;
completed = 0;
workload = 0;
IJ.showProgress(1.0);
} /* end resetProgressBar */
/*********************************************************************
* Perform <code>stride</code> operations at once.
*
* @param stride Amount of work that is skipped.
********************************************************************/
public static synchronized void skipProgressBar (
final int stride
) {
completed += stride - 1;
stepProgressBar();
} /* end skipProgressBar */
/*********************************************************************
* Perform <code>1</code> operation unit.
********************************************************************/
public static synchronized void stepProgressBar (
) {
final long timeStamp = System.currentTimeMillis();
completed = completed + 1;
if ((TIME_QUANTUM <= (timeStamp - lastTime)) | (completed == workload)) {
lastTime = timeStamp;
IJ.showProgress((double)completed / (double)workload);
}
} /* end stepProgressBar */
/*********************************************************************
* Acknowledge that <code>batch</code> work has been performed.
*
* @param batch Completed amount of work.
********************************************************************/
public static synchronized void workloadDone (
final int batch
) {
workload -= batch;
completed -= batch;
} /* end workloadDone */
} /* end class unwarpJProgressBar */
/*====================================================================
| unwarpJTransformation
\===================================================================*/
/*------------------------------------------------------------------*/
class unwarpJTransformation
{ /* begin class unwarpJTransformation */
/*....................................................................
Private variables
....................................................................*/
private final double FLT_EPSILON = (double)Float.intBitsToFloat((int)0x33FFFFFF);
private final boolean PYRAMID = true;
private final boolean ORIGINAL = false;
private final int transformationSplineDegree=3;
// Some useful references
private ImagePlus output_ip;
private unwarpJDialog dialog;
// Images
private ImagePlus sourceImp;
private ImagePlus targetImp;
private unwarpJImageModel source;
private unwarpJImageModel target;
// Landmarks
private unwarpJPointHandler sourcePh;
private unwarpJPointHandler targetPh;
// Masks for the images
private unwarpJMask sourceMsk;
private unwarpJMask targetMsk;
// Image size
private int sourceHeight;
private int sourceWidth;
private int targetHeight;
private int targetWidth;
private int targetCurrentHeight;
private int targetCurrentWidth;
private double factorHeight;
private double factorWidth;
// Transformation parameters
private int min_scale_deformation;
private int max_scale_deformation;
private int min_scale_image;
private int outputLevel;
private boolean showMarquardtOptim;
private double divWeight;
private double curlWeight;
private double landmarkWeight;
private double imageWeight;
private double stopThreshold;
private int accurate_mode;
private boolean saveTransf;
private String fn_tnf;
// Transformation estimate
private int intervals;
private double [][]cx;
private double [][]cy;
private double [][]transformation_x;
private double [][]transformation_y;
private unwarpJImageModel swx;
private unwarpJImageModel swy;
// Regularization temporary variables
private double [][]P11;
private double [][]P22;
private double [][]P12;
public ImageStack is = null;
/*....................................................................
Public methods
....................................................................*/
/*------------------------------------------------------------------*/
public void doRegistration (
) {
// This function can only be applied with splines of an odd order
// Bring into consideration the image/coefficients at the smallest scale
source.popFromPyramid();
target.popFromPyramid();
targetCurrentHeight = target.getCurrentHeight();
targetCurrentWidth = target.getCurrentWidth();
factorHeight = target.getFactorHeight();
factorWidth = target.getFactorWidth();
// Ask memory for the transformation coefficients
intervals=(int)Math.pow(2,min_scale_deformation);
cx = new double[intervals+3][intervals+3];
cy = new double[intervals+3][intervals+3];
// Build matrices for computing the regularization
buildRegularizationTemporary(intervals);
// Ask for memory for the residues
final int K;
if (targetPh!=null) K=targetPh.getPoints().size();
else K=0;
double [] dx=new double[K];
double [] dy=new double[K];
computeInitialResidues(dx,dy);
// Compute the affine transformation from the target to the source coordinates
// Notice that this matrix is independent of the scale, but the residues are not
double[][] affineMatrix = null;
if (landmarkWeight==0) affineMatrix=computeAffineMatrix();
else {
affineMatrix=new double[2][3];
affineMatrix[0][0]=affineMatrix[1][1]=1;
affineMatrix[0][1]=affineMatrix[0][2]=0;
affineMatrix[1][0]=affineMatrix[1][2]=0;
}
// Incorporate the affine transformation into the spline coefficient
for (int i= 0; i<intervals + 3; i++) {
final double v = (double)((i - 1) * (targetCurrentHeight - 1)) / (double)intervals;
final double xv = affineMatrix[0][2] + affineMatrix[0][1] * v;
final double yv = affineMatrix[1][2] + affineMatrix[1][1] * v;
for (int j = 0; j < intervals + 3; j++) {
final double u = (double)((j - 1) * (targetCurrentWidth - 1)) / (double)intervals;
cx[i][j] = xv + affineMatrix[0][0] * u;
cy[i][j] = yv + affineMatrix[1][0] * u;
}
}
// Now refine with the different scales
int state; // state=-1 --> Finish
// state= 0 --> Increase deformation detail
// state= 1 --> Increase image detail
// state= 2 --> Do nothing until the finest image scale
if (min_scale_deformation==max_scale_deformation) state=1;
else state=0;
int s=min_scale_deformation;
int step=0;
computeTotalWorkload();
while (state!=-1) {
int currentDepth=target.getCurrentDepth();
// Update the deformation coefficients only in states 0 and 1
if (state==0 || state==1) {
// Update the deformation coefficients with the error of the landmarks
// The following conditional is now useless but it is there to allow
// easy changes like applying the landmarks only in the coarsest deformation
if (s>=min_scale_deformation) {
// Number of intervals at this scale and ask for memory
intervals=(int) Math.pow(2,s);
final double[][] newcx = new double[intervals+3][intervals+3];
final double[][] newcy = new double[intervals+3][intervals+3];
// Compute the residues before correcting at this scale
computeScaleResidues(intervals, cx, cy, dx, dy);
// Compute the coefficients at this scale
boolean underconstrained=true;
if (divWeight==0 && curlWeight==0)
underconstrained=
computeCoefficientsScale(intervals, dx, dy, newcx, newcy);
else
underconstrained=
computeCoefficientsScaleWithRegularization(
intervals, dx, dy, newcx, newcy);
// Incorporate information from the previous scale
if (!underconstrained || (step==0 && landmarkWeight!=0)) {
for (int i=0; i<intervals+3; i++)
for (int j=0; j<intervals+3; j++) {
cx[i][j]+=newcx[i][j];
cy[i][j]+=newcy[i][j];
}
}
}
// Optimize deformation coefficients
if (imageWeight!=0)
optimizeCoeffs(intervals,stopThreshold,cx,cy);
}
// Prepare for next iteration
step++;
switch (state) {
case 0:
// Finer details in the deformation
if (s<max_scale_deformation) {
cx=propagateCoeffsToNextLevel(intervals,cx,1);
cy=propagateCoeffsToNextLevel(intervals,cy,1);
s++;
intervals*=2;
// Prepare matrices for the regularization term
buildRegularizationTemporary(intervals);
if (currentDepth>min_scale_image) state=1;
else state=0;
} else
if (currentDepth>min_scale_image) state=1;
else state=2;
break;
case 1: // Finer details in the image, go on optimizing
case 2: // Finer details in the image, do not optimize
// Compute next state
if (state==1) {
if (s==max_scale_deformation && currentDepth==min_scale_image) state=2;
else if (s==max_scale_deformation) state=1;
else state=0;
} else if (state==2) {
if (currentDepth==0) state=-1;
else state= 2;
}
// Pop another image and prepare the deformation
if (currentDepth!=0) {
double oldTargetCurrentHeight=targetCurrentHeight;
double oldTargetCurrentWidth =targetCurrentWidth;
source.popFromPyramid();
target.popFromPyramid();
targetCurrentHeight = target.getCurrentHeight();
targetCurrentWidth = target.getCurrentWidth();
factorHeight = target.getFactorHeight();
factorWidth = target.getFactorWidth();
// Adapt the transformation to the new image size
double factorY=(targetCurrentHeight-1)/(oldTargetCurrentHeight-1);
double factorX=(targetCurrentWidth -1)/(oldTargetCurrentWidth -1);
for (int i=0; i<intervals+3; i++)
for (int j=0; j<intervals+3; j++) {
cx[i][j]*=factorX;
cy[i][j]*=factorY;
}
// Prepare matrices for the regularization term
buildRegularizationTemporary(intervals);
}
break;
}
// In accurate_mode reduce the stopping threshold for the last iteration
if ((state==0 || state==1) && s==max_scale_deformation &&
currentDepth==min_scale_image+1 && accurate_mode==1)
stopThreshold/=10;
}
// Show results
showTransformation(intervals,cx,cy);
if (saveTransf) saveTransformation(intervals,cx,cy);
} /* end doMultiresolutionElasticTransformation */
/*--------------------------------------------------------------------------*/
public double evaluateImageSimilarity() {
int int3=intervals+3;
int halfM=int3*int3;
int M=halfM*2;
double []x = new double [M];
double []grad = new double [M];
for (int i= 0, p=0; i<intervals + 3; i++)
for (int j = 0; j < intervals + 3; j++, p++) {
x[ p]=cx[i][j];
x[halfM+p]=cy[i][j];
}
if (swx==null) {
swx=new unwarpJImageModel(cx);
swy=new unwarpJImageModel(cy);
swx.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
swy.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
}
if (swx.precomputed_getWidth()!=target.getCurrentWidth()) {
swx.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
swy.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
}
return evaluateSimilarity(x,intervals,grad,true,false);
}
/*------------------------------------------------------------------*/
public void getDeformation(
final double [][]transformation_x,
final double [][]transformation_y) {
computeDeformation(intervals,cx,cy,
transformation_x,transformation_y);
}
/*------------------------------------------------------------------*/
public unwarpJTransformation (
final ImagePlus sourceImp,
final ImagePlus targetImp,
final unwarpJImageModel source,
final unwarpJImageModel target,
final unwarpJPointHandler sourcePh,
final unwarpJPointHandler targetPh,
final unwarpJMask sourceMsk,
final unwarpJMask targetMsk,
final int min_scale_deformation,
final int max_scale_deformation,
final int min_scale_image,
final double divWeight,
final double curlWeight,
final double landmarkWeight,
final double imageWeight,
final double stopThreshold,
final int outputLevel,
final boolean showMarquardtOptim,
final int accurate_mode,
final boolean saveTransf,
final String fn_tnf,
final ImagePlus output_ip,
final unwarpJDialog dialog
) {
this.sourceImp = sourceImp;
this.targetImp = targetImp;
this.source = source;
this.target = target;
this.sourcePh = sourcePh;
this.targetPh = targetPh;
this.sourceMsk = sourceMsk;
this.targetMsk = targetMsk;
this.min_scale_deformation = min_scale_deformation;
this.max_scale_deformation = max_scale_deformation;
this.min_scale_image = min_scale_image;
this.divWeight = divWeight;
this.curlWeight = curlWeight;
this.landmarkWeight = landmarkWeight;
this.imageWeight = imageWeight;
this.stopThreshold = stopThreshold;
this.outputLevel = outputLevel;
this.showMarquardtOptim = showMarquardtOptim;
this.accurate_mode = accurate_mode;
this.saveTransf = saveTransf;
this.fn_tnf = fn_tnf;
this.output_ip = output_ip;
this.dialog = dialog;
sourceWidth = source.getWidth();
sourceHeight = source.getHeight();
targetWidth = target.getWidth();
targetHeight = target.getHeight();
} /* end unwarpJTransformation */
/*------------------------------------------------------------------*/
public void transform(double u, double v, double []xyF) {
final double tu = (u * intervals) / (double)(target.getCurrentWidth () - 1) + 1.0F;
final double tv = (v * intervals) / (double)(target.getCurrentHeight() - 1) + 1.0F;
final boolean ORIGINAL=false;
swx.prepareForInterpolation(tu,tv,ORIGINAL); xyF[0]=swx.interpolateI();
swy.prepareForInterpolation(tu,tv,ORIGINAL); xyF[1]=swy.interpolateI();
}
/*....................................................................
Private methods
....................................................................*/
/*------------------------------------------------------------------*/
private void build_Matrix_B(
int intervals, // Intervals in the deformation
int K, // Number of landmarks
double [][]B // System matrix of the landmark interpolation
) {
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
for (int k = 0; k<K; k++) {
final Point targetPoint = (Point)targetVector.elementAt(k);
double x=factorWidth *(double)targetPoint.x;
double y=factorHeight*(double)targetPoint.y;
final double[] bx = xWeight(x, intervals, true);
final double[] by = yWeight(y, intervals, true);
for (int i=0; i<intervals+3; i++)
for (int j=0; j<intervals+3; j++)
B[k][(intervals+3)*i+j] = by[i] * bx[j];
}
}
/*------------------------------------------------------------------*/
private void build_Matrix_Rq1q2(
int intervals,
double weight,
int q1, int q2,
double [][]R
){build_Matrix_Rq1q2q3q4(intervals,weight,q1,q2,q1,q2,R);}
private void build_Matrix_Rq1q2q3q4(
int intervals,
double weight,
int q1, int q2, int q3, int q4,
double [][]R
){
/* Let's define alpha_q as the q-th derivative of a B-Spline
q n
d B (x)
alpha_q(x)= --------------
q
dx
eta_q1q2(x,s1,s2)=integral_0^Xdim alpha_q1(x/h-s1) alpha_q2(x/h-s2)
*/
double [][]etaq1q3=new double[16][16];
int Ydim=target.getCurrentHeight();
int Xdim=target.getCurrentWidth();
build_Matrix_R_geteta(etaq1q3,q1,q3,Xdim,intervals);
double [][]etaq2q4=null;
if (q2!=q1 || q4!=q3 || Ydim!=Xdim) {
etaq2q4=new double[16][16];
build_Matrix_R_geteta(etaq2q4,q2,q4,Ydim,intervals);
} else etaq2q4=etaq1q3;
int M=intervals+1;
int Mp=intervals+3;
for (int l=-1; l<=M; l++)
for (int k=-1; k<=M; k++)
for (int n=-1; n<=M; n++)
for (int m=-1; m<=M; m++) {
int []ip=new int[2];
int []jp=new int[2];
boolean valid_i=build_Matrix_R_getetaindex(l,n,intervals,ip);
boolean valid_j=build_Matrix_R_getetaindex(k,m,intervals,jp);
if (valid_i && valid_j) {
int mn=(n+1)*Mp+(m+1);
int kl=(l+1)*Mp+(k+1);
R[kl][mn]+=weight*etaq1q3[jp[0]][jp[1]]*etaq2q4[ip[0]][ip[1]];
}
}
}
/*------------------------------------------------------------------*/
private double build_Matrix_R_computeIntegral_aa(
double x0,
double xF,
double s1,
double s2,
double h,
int q1,
int q2
) {
// Computes the following integral
//
// xF d^q1 3 x d^q2 3 x
// integral ----- B (--- - s1) ----- B (--- - s2) dx
// x0 dx^q1 h dx^q2 h
// Form the spline coefficients
double [][]C=new double [3][3];
int [][]d=new int [3][3];
double [][]s=new double [3][3];
C[0][0]= 1 ; C[0][1]= 0 ; C[0][2]= 0 ;
C[1][0]= 1 ; C[1][1]=-1 ; C[1][2]= 0 ;
C[2][0]= 1 ; C[2][1]=-2 ; C[2][2]= 1 ;
d[0][0]= 3 ; d[0][1]= 0 ; d[0][2]= 0 ;
d[1][0]= 2 ; d[1][1]= 2 ; d[1][2]= 0 ;
d[2][0]= 1 ; d[2][1]= 1 ; d[2][2]= 1 ;
s[0][0]= 0 ; s[0][1]= 0 ; s[0][2]= 0 ;
s[1][0]=-0.5; s[1][1]= 0.5; s[1][2]= 0 ;
s[2][0]= 1 ; s[2][1]= 0 ; s[2][2]=-1 ;
// Compute the integral
double integral=0;
for (int k=0; k<3; k++) {
double ck=C[q1][k]; if (ck==0) continue;
for (int l=0; l<3; l++) {
double cl=C[q2][l]; if (cl==0) continue;
integral+=ck*cl*build_matrix_R_computeIntegral_BB(
x0,xF,s1+s[q1][k],s2+s[q2][l],h,d[q1][k],d[q2][l]);
}
}
return integral;
}
/*------------------------------------------------------------------*/
private double build_matrix_R_computeIntegral_BB(
double x0,
double xF,
double s1,
double s2,
double h,
int n1,
int n2
) {
// Computes the following integral
//
// xF n1 x n2 x
// integral B (--- - s1) B (--- - s2) dx
// x0 h h
// Change the variable so that the h disappears
// X=x/h
double xFp=xF/h;
double x0p=x0/h;
// Form the spline coefficients
double []c1=new double [n1+2];
double fact_n1=1; for (int k=2; k<=n1; k++) fact_n1*=k;
double sign=1;
for (int k=0; k<=n1+1; k++, sign*=-1)
c1[k]=sign*unwarpJMathTools.nchoosek(n1+1,k)/fact_n1;
double []c2=new double [n2+2];
double fact_n2=1; for (int k=2; k<=n2; k++) fact_n2*=k;
sign=1;
for (int k=0; k<=n2+1; k++, sign*=-1)
c2[k]=sign*unwarpJMathTools.nchoosek(n2+1,k)/fact_n2;
// Compute the integral
double n1_2=(double)((n1+1))/2.0;
double n2_2=(double)((n2+1))/2.0;
double integral=0;
for (int k=0; k<=n1+1; k++)
for (int l=0; l<=n2+1; l++) {
integral+=
c1[k]*c2[l]*build_matrix_R_computeIntegral_xx(
x0p,xFp,s1+k-n1_2,s2+l-n2_2,n1,n2);
}
return integral*h;
}
/*------------------------------------------------------------------*/
private double build_matrix_R_computeIntegral_xx(
double x0,
double xF,
double s1,
double s2,
int q1,
int q2
){
// Computation of the integral
// xF q1 q2
// integral (x-s1) (x-s2) dx
// x0 + +
// Change of variable so that s1 is 0
// X=x-s1 => x-s2=X-(s2-s1)
double s2p=s2-s1;
double xFp=xF-s1;
double x0p=x0-s1;
// Now integrate
if (xFp<0) return 0;
// Move x0 to the first point where both integrands
// are distinct from 0
x0p=Math.max(x0p,Math.max(s2p,0));
if (x0p>xFp) return 0;
// There is something to integrate
// Evaluate the primitive at xF and x0
double IxFp=0;
double Ix0p=0;
for (int k=0; k<=q2; k++) {
double aux=unwarpJMathTools.nchoosek(q2,k)/(q1+k+1)*
Math.pow(-s2p,q2-k);
IxFp+=Math.pow(xFp,q1+k+1)*aux;
Ix0p+=Math.pow(x0p,q1+k+1)*aux;
}
return IxFp-Ix0p;
}
/*------------------------------------------------------------------*/
private void build_Matrix_R_geteta(
double [][]etaq1q2,
int q1,
int q2,
int dim,
int intervals
) {
boolean [][]done=new boolean[16][16];
// Clear
for (int i=0; i<16; i++)
for (int j=0; j<16; j++) {
etaq1q2[i][j]=0;
done[i][j]=false;
}
// Compute each integral needed
int M=intervals+1;
double h=(double)dim/intervals;
for (int ki1=-1; ki1<=M; ki1++)
for (int ki2=-1; ki2<=M; ki2++) {
int []ip=new int[2];
boolean valid_i=build_Matrix_R_getetaindex(ki1,ki2,intervals,ip);
if (valid_i && !done[ip[0]][ip[1]]) {
etaq1q2[ip[0]][ip[1]]=
build_Matrix_R_computeIntegral_aa(0,dim,ki1,ki2,h,q1,q2);
done[ip[0]][ip[1]]=true;
}
}
}
/*------------------------------------------------------------------*/
private boolean build_Matrix_R_getetaindex(
int ki1,
int ki2,
int intervals,
int []ip
) {
ip[0]=0;
ip[1]=0;
// Determine the clipped inner limits of the intersection
int kir=Math.min(intervals,Math.min(ki1,ki2)+2);
int kil=Math.max(0 ,Math.max(ki1,ki2)-2);
if (kil>=kir) return false;
// Determine which are the pieces of the
// function that lie in the intersection
int two_i=1;
double ki;
for (int i=0; i<=3; i++, two_i*=2) {
// First function
ki=ki1+i-1.5; // Middle sample of the piece i
if (kil<=ki && ki<=kir) ip[0]+=two_i;
// Second function
ki=ki2+i-1.5; // Middle sample of the piece i
if (kil<=ki && ki<=kir) ip[1]+=two_i;
}
ip[0]--;
ip[1]--;
return true;
}
/*------------------------------------------------------------------*/
private void buildRegularizationTemporary(int intervals) {
// M is the number of spline coefficients per row
int M=intervals+3;
int M2=M*M;
// P11
P11=new double[M2][M2];
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) P11[i][j]=0.0;
build_Matrix_Rq1q2(intervals, divWeight , 2, 0, P11);
build_Matrix_Rq1q2(intervals, divWeight+curlWeight, 1, 1, P11);
build_Matrix_Rq1q2(intervals, curlWeight, 0, 2, P11);
// P22
P22=new double[M2][M2];
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) P22[i][j]=0.0;
build_Matrix_Rq1q2(intervals, divWeight , 0, 2, P22);
build_Matrix_Rq1q2(intervals, divWeight+curlWeight, 1, 1, P22);
build_Matrix_Rq1q2(intervals, curlWeight, 2, 0, P22);
// P12
P12=new double[M2][M2];
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) P12[i][j]=0.0;
build_Matrix_Rq1q2q3q4(intervals, 2*divWeight , 2, 0, 1, 1, P12);
build_Matrix_Rq1q2q3q4(intervals, 2*divWeight , 1, 1, 0, 2, P12);
build_Matrix_Rq1q2q3q4(intervals,-2*curlWeight, 0, 2, 1, 1, P12);
build_Matrix_Rq1q2q3q4(intervals,-2*curlWeight, 1, 1, 2, 0, P12);
}
/*------------------------------------------------------------------*/
private double[][] computeAffineMatrix (
) {
boolean adjust_size=false;
final double[][] D = new double[3][3];
final double[][] H = new double[3][3];
final double[][] U = new double[3][3];
final double[][] V = new double[3][3];
final double[][] X = new double[2][3];
final double[] W = new double[3];
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
int removeLastPoint=0;
if (false) {
removeLastPoint=sourceMsk.numberOfMaskPoints();
for (int i=0; i<removeLastPoint; i++) {
sourceVector.addElement(sourceMsk.getPoint(i));
targetVector.addElement(targetMsk.getPoint(i));
}
}
int n = targetVector.size();
switch (n) {
case 0:
for (int i = 0; (i < 2); i++)
for (int j = 0; (j < 3); j++) X[i][j]=0.0;
if (adjust_size) {
// Make both images of the same size
X[0][0]=(double)source.getCurrentWidth ()/target.getCurrentWidth ();
X[1][1]=(double)source.getCurrentHeight()/target.getCurrentHeight();
} else {
// Make both images to be centered
X[0][0]=X[1][1]=1;
X[0][2]=((double)source.getCurrentWidth ()-target.getCurrentWidth ())/2;
X[1][2]=((double)source.getCurrentHeight()-target.getCurrentHeight())/2;
}
break;
case 1:
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 2); j++) {
X[i][j] = (i == j) ? (1.0F) : (0.0F);
}
}
X[0][2] = factorWidth*(double)(((Point)sourceVector.firstElement()).x
- ((Point)targetVector.firstElement()).x);
X[1][2] = factorHeight*(double)(((Point)sourceVector.firstElement()).y
- ((Point)targetVector.firstElement()).y);
break;
case 2:
final double x0 = factorWidth *((Point)sourceVector.elementAt(0)).x;
final double y0 = factorHeight*((Point)sourceVector.elementAt(0)).y;
final double x1 = factorWidth *((Point)sourceVector.elementAt(1)).x;
final double y1 = factorHeight*((Point)sourceVector.elementAt(1)).y;
final double u0 = factorWidth *((Point)targetVector.elementAt(0)).x;
final double v0 = factorHeight*((Point)targetVector.elementAt(0)).y;
final double u1 = factorWidth *((Point)targetVector.elementAt(1)).x;
final double v1 = factorHeight*((Point)targetVector.elementAt(1)).y;
sourceVector.addElement(new Point((int)(x1 + y0 - y1), (int)(x1 + y1 - x0)));
targetVector.addElement(new Point((int)(u1 + v0 - v1), (int)(u1 + v1 - u0)));
removeLastPoint=1;
n = 3;
default:
for (int i = 0; (i < 3); i++) {
for (int j = 0; (j < 3); j++) {
H[i][j] = 0.0F;
}
}
for (int k = 0; (k < n); k++) {
final Point sourcePoint = (Point)sourceVector.elementAt(k);
final Point targetPoint = (Point)targetVector.elementAt(k);
final double sx=factorWidth * (double)sourcePoint.x;
final double sy=factorHeight* (double)sourcePoint.y;
final double tx=factorWidth * (double)targetPoint.x;
final double ty=factorHeight* (double)targetPoint.y;
H[0][0] += tx * sx;
H[0][1] += tx * sy;
H[0][2] += tx;
H[1][0] += ty * sx;
H[1][1] += ty * sy;
H[1][2] += ty;
H[2][0] += sx;
H[2][1] += sy;
H[2][2] += 1.0F;
D[0][0] += sx * sx;
D[0][1] += sx * sy;
D[0][2] += sx;
D[1][0] += sy * sx;
D[1][1] += sy * sy;
D[1][2] += sy;
D[2][0] += sx;
D[2][1] += sy;
D[2][2] += 1.0F;
}
unwarpJMathTools.singularValueDecomposition(H, W, V);
if ((Math.abs(W[0]) < FLT_EPSILON) || (Math.abs(W[1]) < FLT_EPSILON)
|| (Math.abs(W[2]) < FLT_EPSILON)) {
return(computeRotationMatrix());
}
for (int i = 0; (i < 3); i++) {
for (int j = 0; (j < 3); j++) {
V[i][j] /= W[j];
}
}
for (int i = 0; (i < 3); i++) {
for (int j = 0; (j < 3); j++) {
U[i][j] = 0.0F;
for (int k = 0; (k < 3); k++) {
U[i][j] += D[i][k] * V[k][j];
}
}
}
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 3); j++) {
X[i][j] = 0.0F;
for (int k = 0; (k < 3); k++) {
X[i][j] += U[i][k] * H[j][k];
}
}
}
break;
}
if (removeLastPoint!=0) {
for (int i=1; i<=removeLastPoint; i++) {
sourcePh.getPoints().removeElementAt(n-i);
targetPh.getPoints().removeElementAt(n-i);
}
}
return(X);
} /* end computeAffineMatrix */
/*------------------------------------------------------------------*/
private void computeAffineResidues(
final double[][] affineMatrix, // Input
final double[] dx, // output, difference in x for each landmark
final double[] dy // output, difference in y for each landmark
// The output vectors should be already resized
) {
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
final int K = targetPh.getPoints().size();
for (int k=0; k<K; k++) {
final Point sourcePoint = (Point)sourceVector.elementAt(k);
final Point targetPoint = (Point)targetVector.elementAt(k);
double u = factorWidth *(double)targetPoint.x;
double v = factorHeight*(double)targetPoint.y;
final double x = affineMatrix[0][2]
+ affineMatrix[0][0] * u + affineMatrix[0][1] * v;
final double y = affineMatrix[1][2]
+ affineMatrix[1][0] * u + affineMatrix[1][1] * v;
dx[k] = factorWidth*(double)sourcePoint.x - x;
dy[k] = factorHeight*(double)sourcePoint.y - y;
}
}
/*------------------------------------------------------------------*/
private boolean computeCoefficientsScale(
final int intervals, // input, number of intervals at this scale
final double []dx, // input, x residue so far
final double []dy, // input, y residue so far
final double [][]cx, // output, x coefficients for splines
final double [][]cy // output, y coefficients for splines
) {
int K=0;
if (targetPh!=null) K=targetPh.getPoints().size();
boolean underconstrained=false;
if (0<K) {
// Form the equation system Bc=d
final double[][] B = new double[K][(intervals + 3) * (intervals + 3)];
build_Matrix_B(intervals, K, B);
// "Invert" the matrix B
int Nunk=(intervals+3)*(intervals+3);
double [][] iB=new double[Nunk][K];
underconstrained=unwarpJMathTools.invertMatrixSVD(K,Nunk,B,iB);
// Now multiply iB times dx and dy respectively
int ij=0;
for (int i = 0; i<intervals+3; i++)
for (int j = 0; j<intervals+3; j++) {
cx[i][j] = cy[i][j] = 0.0F;
for (int k = 0; k<K; k++) {
cx[i][j] += iB[ij][k] * dx[k];
cy[i][j] += iB[ij][k] * dy[k];
}
ij++;
}
}
return underconstrained;
}
/*------------------------------------------------------------------*/
private boolean computeCoefficientsScaleWithRegularization(
final int intervals, // input, number of intervals at this scale
final double []dx, // input, x residue so far
final double []dy, // input, y residue so far
final double [][]cx, // output, x coefficients for splines
final double [][]cy // output, y coefficients for splines
) {
boolean underconstrained=true;
int K = 0;
if (targetPh!=null) K=targetPh.getPoints().size();
if (0<K) {
// M is the number of spline coefficients per row
int M=intervals+3;
int M2=M*M;
// Create A and b for the system Ac=b
final double[][] A=new double[2*M2][2*M2];
final double[] b=new double[2*M2];
for (int i=0; i<2*M2; i++) {
b[i]=0.0;
for (int j=0; j<2*M2; j++) A[i][j]=0.0;
}
// Get the matrix related to the landmarks
final double[][] B = new double[K][M2];
build_Matrix_B(intervals, K, B);
// Fill the part of the equation system related to the landmarks
// Compute 2 * B^t * B
for (int i=0; i<M2; i++) {
for (int j=i; j<M2; j++) {
double bitbj=0; // bi^t * bj, i.e., column i x column j
for (int l=0; l<K; l++)
bitbj+=B[l][i]*B[l][j];
bitbj*=2;
int ij=i*M2+j;
A[M2+i][M2+j]=A[M2+j][M2+i]=A[i][j]=A[j][i]=bitbj;
}
}
// Compute 2 * B^t * [dx dy]
for (int i=0; i<M2; i++) {
double bitdx=0;
double bitdy=0;
for (int l=0; l<K; l++) {
bitdx+=B[l][i]*dx[l];
bitdy+=B[l][i]*dy[l];
}
bitdx*=2;
bitdy*=2;
b[ i]=bitdx;
b[M2+i]=bitdy;
}
// Get the matrices associated to the regularization
// Copy P11 symmetrized to the equation system
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) {
double aux=P11[i][j];
A[i][j]+=aux;
A[j][i]+=aux;
}
// Copy P22 symmetrized to the equation system
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) {
double aux=P22[i][j];
A[M2+i][M2+j]+=aux;
A[M2+j][M2+i]+=aux;
}
// Copy P12 and P12^t to their respective places
for (int i=0; i<M2; i++)
for (int j=0; j<M2; j++) {
A[ i][M2+j]=P12[i][j]; // P12
A[M2+i][ j]=P12[j][i]; // P12^t
}
// Now solve the system
// Invert the matrix A
double [][] iA=new double[2*M2][2*M2];
underconstrained=unwarpJMathTools.invertMatrixSVD(2*M2,2*M2,A,iA);
// Now multiply iB times b and distribute in cx and cy
int ij=0;
for (int i = 0; i<intervals+3; i++)
for (int j = 0; j<intervals+3; j++) {
cx[i][j] = cy[i][j] = 0.0F;
for (int l = 0; l<2*M2; l++) {
cx[i][j] += iA[ ij][l] * b[l];
cy[i][j] += iA[M2+ij][l] * b[l];
}
ij++;
}
}
return underconstrained;
}
/*------------------------------------------------------------------*/
private void computeInitialResidues(
final double[] dx, // output, difference in x for each landmark
final double[] dy // output, difference in y for each landmark
// The output vectors should be already resized
) {
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
int K = 0;
if (targetPh!=null) targetPh.getPoints().size();
for (int k=0; k<K; k++) {
final Point sourcePoint = (Point)sourceVector.elementAt(k);
final Point targetPoint = (Point)targetVector.elementAt(k);
dx[k] = factorWidth *(sourcePoint.x - targetPoint.x);
dy[k] = factorHeight*(sourcePoint.y - targetPoint.y);
}
}
/*------------------------------------------------------------------*/
private void computeDeformation(
final int intervals,
final double [][]cx,
final double [][]cy,
final double [][]transformation_x,
final double [][]transformation_y) {
// Set these coefficients to an interpolator
unwarpJImageModel swx = new unwarpJImageModel(cx);
unwarpJImageModel swy = new unwarpJImageModel(cy);
// Compute the transformation mapping
for (int v=0; v<targetCurrentHeight; v++) {
final double tv = (double)(v * intervals) / (double)(targetCurrentHeight - 1) + 1.0F;
for (int u = 0; u<targetCurrentWidth; u++) {
final double tu = (double)(u * intervals) / (double)(targetCurrentWidth - 1) + 1.0F;
swx.prepareForInterpolation(tu, tv, ORIGINAL); transformation_x[v][u] = swx.interpolateI();
swy.prepareForInterpolation(tu, tv, ORIGINAL); transformation_y[v][u] = swy.interpolateI();
}
}
}
/*------------------------------------------------------------------*/
private double[][] computeRotationMatrix (
) {
final double[][] X = new double[2][3];
final double[][] H = new double[2][2];
final double[][] V = new double[2][2];
final double[] W = new double[2];
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
final int n = targetVector.size();
switch (n) {
case 0:
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 3); j++) {
X[i][j] = (i == j) ? (1.0F) : (0.0F);
}
}
break;
case 1:
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 2); j++) {
X[i][j] = (i == j) ? (1.0F) : (0.0F);
}
}
X[0][2] = factorWidth *(double)(((Point)sourceVector.firstElement()).x
- ((Point)targetVector.firstElement()).x);
X[1][2] = factorHeight*(double)(((Point)sourceVector.firstElement()).y
- ((Point)targetVector.firstElement()).y);
break;
default:
double xTargetAverage = 0.0F;
double yTargetAverage = 0.0F;
for (int i = 0; (i < n); i++) {
final Point p = (Point)targetVector.elementAt(i);
xTargetAverage += factorWidth *(double)p.x;
yTargetAverage += factorHeight*(double)p.y;
}
xTargetAverage /= (double)n;
yTargetAverage /= (double)n;
final double[] xCenteredTarget = new double[n];
final double[] yCenteredTarget = new double[n];
for (int i = 0; (i < n); i++) {
final Point p = (Point)targetVector.elementAt(i);
xCenteredTarget[i] = factorWidth *(double)p.x - xTargetAverage;
yCenteredTarget[i] = factorHeight*(double)p.y - yTargetAverage;
}
double xSourceAverage = 0.0F;
double ySourceAverage = 0.0F;
for (int i = 0; (i < n); i++) {
final Point p = (Point)sourceVector.elementAt(i);
xSourceAverage += factorWidth *(double)p.x;
ySourceAverage += factorHeight*(double)p.y;
}
xSourceAverage /= (double)n;
ySourceAverage /= (double)n;
final double[] xCenteredSource = new double[n];
final double[] yCenteredSource = new double[n];
for (int i = 0; (i < n); i++) {
final Point p = (Point)sourceVector.elementAt(i);
xCenteredSource[i] = factorWidth *(double)p.x - xSourceAverage;
yCenteredSource[i] = factorHeight*(double)p.y - ySourceAverage;
}
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 2); j++) {
H[i][j] = 0.0F;
}
}
for (int k = 0; (k < n); k++) {
H[0][0] += xCenteredTarget[k] * xCenteredSource[k];
H[0][1] += xCenteredTarget[k] * yCenteredSource[k];
H[1][0] += yCenteredTarget[k] * xCenteredSource[k];
H[1][1] += yCenteredTarget[k] * yCenteredSource[k];
}
// COSS: Watch out that this H is the transpose of the one
// defined in the text. That is why X=V*U^t is the inverse of
// of the rotation matrix.
unwarpJMathTools.singularValueDecomposition(H, W, V);
if (((H[0][0] * H[1][1] - H[0][1] * H[1][0])
* (V[0][0] * V[1][1] - V[0][1] * V[1][0])) < 0.0F) {
if (W[0] < W[1]) {
V[0][0] *= -1.0F;
V[1][0] *= -1.0F;
}
else {
V[0][1] *= -1.0F;
V[1][1] *= -1.0F;
}
}
for (int i = 0; (i < 2); i++) {
for (int j = 0; (j < 2); j++) {
X[i][j] = 0.0F;
for (int k = 0; (k < 2); k++) {
X[i][j] += V[i][k] * H[j][k];
}
}
}
X[0][2] = xSourceAverage - X[0][0] * xTargetAverage - X[0][1] * yTargetAverage;
X[1][2] = ySourceAverage - X[1][0] * xTargetAverage - X[1][1] * yTargetAverage;
break;
}
return(X);
} /* end computeRotationMatrix */
/*------------------------------------------------------------------*/
private void computeScaleResidues(
int intervals, // input, number of intevals
final double [][]cx, // Input, spline coefficients
final double [][]cy,
final double []dx, // Input/Output. At the input it has the
final double []dy // residue so far, at the output this
// residue is modified to account for
// the model at the new scale
) {
// Set these coefficients to an interpolator
unwarpJImageModel swx = new unwarpJImageModel(cx);
unwarpJImageModel swy = new unwarpJImageModel(cy);
// Get the list of landmarks
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
final int K = targetVector.size();
for (int k=0; k<K; k++) {
// Get the landmark coordinate in the target image
final Point sourcePoint = (Point)sourceVector.elementAt(k);
final Point targetPoint = (Point)targetVector.elementAt(k);
double u = factorWidth *(double)targetPoint.x;
double v = factorHeight*(double)targetPoint.y;
// Express it in "spline" units
double tu = (double)(u * intervals) / (double)(targetCurrentWidth - 1) + 1.0F;
double tv = (double)(v * intervals) / (double)(targetCurrentHeight - 1) + 1.0F;
// Transform this coordinate to the source image
swx.prepareForInterpolation(tu, tv, false); double x=swx.interpolateI();
swy.prepareForInterpolation(tu, tv, false); double y=swy.interpolateI();
// Substract the result from the residual
dx[k]=factorWidth *(double)sourcePoint.x - x;
dy[k]=factorHeight*(double)sourcePoint.y - y;
}
}
/*--------------------------------------------------------------------------*/
private void computeTotalWorkload() {
// This code is an excerpt from doRegistration() to compute the exact
// number of steps
// Now refine with the different scales
int state; // state=-1 --> Finish
// state= 0 --> Increase deformation detail
// state= 1 --> Increase image detail
// state= 2 --> Do nothing until the finest image scale
if (min_scale_deformation==max_scale_deformation) state=1;
else state=0;
int s=min_scale_deformation;
int currentDepth=target.getCurrentDepth();
int workload=0;
while (state!=-1) {
// Update the deformation coefficients only in states 0 and 1
if (state==0 || state==1) {
// Optimize deformation coefficients
if (imageWeight!=0)
workload+=300*(currentDepth+1);
}
// Prepare for next iteration
switch (state) {
case 0:
// Finer details in the deformation
if (s<max_scale_deformation) {
s++;
if (currentDepth>min_scale_image) state=1;
else state=0;
} else
if (currentDepth>min_scale_image) state=1;
else state=2;
break;
case 1: // Finer details in the image, go on optimizing
case 2: // Finer details in the image, do not optimize
// Compute next state
if (state==1) {
if (s==max_scale_deformation && currentDepth==min_scale_image) state=2;
else if (s==max_scale_deformation) state=1;
else state=0;
} else if (state==2) {
if (currentDepth==0) state=-1;
else state= 2;
}
// Pop another image and prepare the deformation
if (currentDepth!=0) currentDepth--;
break;
}
}
unwarpJProgressBar.resetProgressBar();
unwarpJProgressBar.addWorkload(workload);
}
/*--------------------------------------------------------------------------*/
private double evaluateSimilarity(
final double []c, // Input: Deformation coefficients
final int intervals, // Input: Number of intervals for the deformation
double []grad, // Output: Gradient of the similarity
// Output: the similarity is returned
final boolean only_image, // Input: if true, only the image term is considered
// and not the regularization
final boolean show_error // Input: if true, an image is shown with the error
) {
int cYdim=intervals+3;
int cXdim=cYdim;
int Nk=cYdim*cXdim;
int twiceNk=2*Nk;
double []vgradreg =new double[grad.length];
double []vgradland=new double[grad.length];
// Set the transformation coefficients to the interpolator
swx.setCoefficients(c,cYdim,cXdim,0);
swy.setCoefficients(c,cYdim,cXdim,Nk);
// Initialize gradient
for (int k=0; k<twiceNk; k++) vgradreg[k]=vgradland[k]=grad[k]=0.0F;
// Estimate the similarity and gradient between both images
double imageSimilarity=0.0;
int Ydim=target.getCurrentHeight();
int Xdim=target.getCurrentWidth();
// Prepare to show
double [][]error_image=null;
double [][]div_error_image=null;
double [][]curl_error_image=null;
double [][]laplacian_error_image=null;
double [][]jacobian_error_image=null;
if (show_error) {
error_image=new double[Ydim][Xdim];
div_error_image=new double[Ydim][Xdim];
curl_error_image=new double[Ydim][Xdim];
laplacian_error_image=new double[Ydim][Xdim];
jacobian_error_image=new double[Ydim][Xdim];
for (int v=0; v<Ydim; v++)
for (int u=0; u<Xdim; u++)
error_image[v][u]=div_error_image[v][u]=curl_error_image[v][u]=
laplacian_error_image[v][u]=jacobian_error_image[v][u]=-1.0;
}
// Loop over all points in the source image
int n=0;
if (imageWeight!=0 || show_error) {
double []xD2=new double[3]; // Some space for the second derivatives
double []yD2=new double[3]; // of the transformation
double []xD =new double[2]; // Some space for the second derivatives
double []yD =new double[2]; // of the transformation
double []I1D=new double[2]; // Space for the first derivatives of I1
double hx=(Xdim-1)/intervals; // Scale in the X axis
double hy=(Ydim-1)/intervals; // Scale in the Y axis
double []targetCurrentImage=target.getCurrentImage();
int uv=0;
for (int v=0; v<Ydim; v++) {
for (int u=0; u<Xdim; u++, uv++) {
// Compute image term .....................................................
// Check if this point is in the target mask
if (targetMsk.getValue(u/factorWidth,v/factorHeight)) {
// Compute value in the source image
double I2=targetCurrentImage[uv];
// Compute the position of this point in the target
double x=swx.precomputed_interpolateI(u,v);
double y=swy.precomputed_interpolateI(u,v);
// Check if this point is in the source mask
if (sourceMsk.getValue(x/factorWidth,y/factorHeight)) {
// Compute the value of the target at that point
source.prepareForInterpolation(x,y,PYRAMID); double I1=source.interpolateI();
source.interpolateD(I1D); double I1dx=I1D[0], I1dy=I1D[1];
double error=I2-I1;
double error2=error*error;
if (show_error) error_image[v][u]=error;
imageSimilarity+=error2;
// Compute the derivative with respect to all the c coefficients
// Cost of the derivatives = 16*(3 mults + 2 sums)
// Current cost= 359 mults + 346 sums
for (int l=0; l<4; l++)
for (int m=0; m<4; m++) {
if (swx.prec_yIndex[v][l]==-1 || swx.prec_xIndex[u][m]==-1) continue;
double weightI=swx.precomputed_getWeightI(l,m,u,v);
int k=swx.prec_yIndex[v][l]*cYdim+swx.prec_xIndex[u][m];
// Compute partial result
// There's also a multiplication by 2 that I will
// do later
double aux=-error*weightI;
// Derivative related to X deformation
grad[k] +=aux*I1dx;
// Derivative related to Y deformation
grad[k+Nk]+=aux*I1dy;
}
n++; // Another point has been successfully evaluated
}
}
// Show regularization images ...........................................
if (show_error) {
double gradcurlx=0.0, gradcurly=0.0;
double graddivx =0.0, graddivy =0.0;
double xdx =0.0, xdy =0.0,
ydx =0.0, ydy =0.0,
xdxdy=0.0, xdxdx=0.0, xdydy=0.0,
ydxdy=0.0, ydxdx=0.0, ydydy=0.0;
// Compute the first derivative terms
swx.precomputed_interpolateD(xD,u,v); xdx=xD[0]/hx; xdy=xD[1]/hy;
swy.precomputed_interpolateD(yD,u,v); ydx=yD[0]/hx; ydy=yD[1]/hy;
// Compute the second derivative terms
swx.precomputed_interpolateD2(xD2,u,v);
xdxdy=xD2[0]; xdxdx=xD2[1]; xdydy=xD2[2];
swy.precomputed_interpolateD2(yD2,u,v);
ydxdy=yD2[0]; ydxdx=yD2[1]; ydydy=yD2[2];
// Error in the divergence
graddivx=xdxdx+ydxdy;
graddivy=xdxdy+ydydy;
double graddiv=graddivx*graddivx+graddivy*graddivy;
double errorgraddiv=divWeight*graddiv;
if (divWeight!=0) div_error_image [v][u]=errorgraddiv;
else div_error_image [v][u]=graddiv;
// Compute error in the curl
gradcurlx=-xdxdy+ydxdx;
gradcurly=-xdydy+ydxdy;
double gradcurl=gradcurlx*gradcurlx+gradcurly*gradcurly;
double errorgradcurl=curlWeight*gradcurl;
if (curlWeight!=0) curl_error_image[v][u]=errorgradcurl;
else curl_error_image[v][u]=gradcurl;
// Compute Laplacian error
laplacian_error_image[v][u] =xdxdx*xdxdx;
laplacian_error_image[v][u]+=xdxdy*xdxdy;
laplacian_error_image[v][u]+=xdydy*xdydy;
laplacian_error_image[v][u]+=ydxdx*ydxdx;
laplacian_error_image[v][u]+=ydxdy*ydxdy;
laplacian_error_image[v][u]+=ydydy*ydydy;
// Compute jacobian error
jacobian_error_image[v][u] =xdx*ydy-xdy*ydx;
}
}
}
}
// Average the image related terms
if (n!=0) {
imageSimilarity*=imageWeight/n;
double aux=imageWeight*2.0/n; // This is the 2 coming from the
// derivative that I would do later
for (int k=0; k<twiceNk; k++) grad[k]*=aux;
} else
if (imageWeight==0) imageSimilarity=0;
else imageSimilarity=1/FLT_EPSILON;
// Compute regularization term ..............................................
double regularization=0.0;
if (!only_image) {
for (int i=0; i<Nk; i++)
for (int j=0; j<Nk; j++) {
regularization+=c[ i]*P11[i][j]*c[ j]+// c1^t P11 c1
c[Nk+i]*P22[i][j]*c[Nk+j]+// c2^t P22 c2
c[ i]*P12[i][j]*c[Nk+j];// c1^t P12 c2
vgradreg[ i]+=2*P11[i][j]*c[j]; // 2 P11 c1
vgradreg[Nk+i]+=2*P22[i][j]*c[Nk+j]; // 2 P22 c2
vgradreg[ i]+= P12[i][j]*c[Nk+j]; // P12 c2
vgradreg[Nk+i]+= P12[j][i]*c[ j]; // P12^t c1
}
regularization*=1.0/(Ydim*Xdim);
for (int k=0; k<twiceNk; k++) vgradreg [k]*=1.0/(Ydim*Xdim);
}
// Compute landmark error and derivative ...............................
// Get the list of landmarks
double landmarkError=0.0;
int K = 0;
if (targetPh!=null) K=targetPh.getPoints().size();
if (landmarkWeight!=0) {
Vector sourceVector=null;
if (sourcePh!=null) sourceVector=sourcePh.getPoints();
else sourceVector=new Vector();
Vector targetVector = null;
if (targetPh!=null) targetVector=targetPh.getPoints();
else targetVector=new Vector();
for (int kp=0; kp<K; kp++) {
// Get the landmark coordinate in the target image
final Point sourcePoint = (Point)sourceVector.elementAt(kp);
final Point targetPoint = (Point)targetVector.elementAt(kp);
double u = factorWidth *(double)targetPoint.x;
double v = factorHeight*(double)targetPoint.y;
// Express it in "spline" units
double tu = (double)(u * intervals) / (double)(targetCurrentWidth - 1) + 1.0F;
double tv = (double)(v * intervals) / (double)(targetCurrentHeight - 1) + 1.0F;
// Transform this coordinate to the source image
swx.prepareForInterpolation(tu, tv, false); double x=swx.interpolateI();
swy.prepareForInterpolation(tu, tv, false); double y=swy.interpolateI();
// Substract the result from the residual
double dx=factorWidth *(double)sourcePoint.x - x;
double dy=factorHeight*(double)sourcePoint.y - y;
// Add to landmark error
landmarkError+=dx*dx+dy*dy;
// Compute the derivative with respect to all the c coefficients
for (int l=0; l<4; l++)
for (int m=0; m<4; m++) {
if (swx.yIndex[l]==-1 || swx.xIndex[m]==-1) continue;
int k=swx.yIndex[l]*cYdim+swx.xIndex[m];
// There's also a multiplication by 2 that I will do later
// Derivative related to X deformation
vgradland[k] -=dx*swx.getWeightI(l,m);
// Derivative related to Y deformation
vgradland[k+Nk]-=dy*swy.getWeightI(l,m);
}
}
}
if (K!=0) {
landmarkError*=landmarkWeight/K;
double aux=2.0*landmarkWeight/K;
// This is the 2 coming from the derivative
// computation that I would do at the end
for (int k=0; k<twiceNk; k++) vgradland[k]*=aux;
}
if (only_image) landmarkError=0;
// Finish computations .............................................................
// Add all gradient terms
for (int k=0; k<twiceNk; k++)
grad[k]+=vgradreg[k]+vgradland[k];
if (show_error) {
unwarpJMiscTools.showImage("Error",error_image);
unwarpJMiscTools.showImage("Divergence Error",div_error_image);
unwarpJMiscTools.showImage("Curl Error",curl_error_image);
unwarpJMiscTools.showImage("Laplacian Error",laplacian_error_image);
unwarpJMiscTools.showImage("Jacobian Error",jacobian_error_image);
}
if (showMarquardtOptim) {
if (imageWeight!=0) IJ.write(" Image error:"+imageSimilarity);
if (landmarkWeight!=0) IJ.write(" Landmark error:"+landmarkError);
if (divWeight!=0 || curlWeight!=0) IJ.write(" Regularization error:"+regularization);
}
return imageSimilarity+landmarkError+regularization;
}
/*--------------------------------------------------------------------------*/
private void Marquardt_it (
double []x,
boolean []optimize,
double []gradient,
double []Hessian,
double lambda
)
{
/* In this function the system (H+lambda*Diag(H))*update=gradient
is solved for update.
H is the hessian of the function f,
gradient is the gradient of the function f,
Diag(H) is a matrix with the diagonal of H.
*/
final double TINY = FLT_EPSILON;
final int M = x.length;
final int Mmax = 35;
final int Mused = Math.min(M,Mmax);
double [][] u = new double [Mused][Mused];
double [][] v = null; //new double [Mused][Mused];
double [] w = null; //new double [Mused];
double [] g = new double [Mused];
double [] update = new double [Mused];
boolean [] optimizep = new boolean [M];
System.arraycopy(optimize,0,optimizep,0,M);
lambda+=1.0F;
if (M>Mmax) {
/* Find the threshold for the most important components */
double []sortedgradient= new double [M];
for (int i=0; i<M; i++) sortedgradient[i]=Math.abs(gradient[i]);
Arrays.sort(sortedgradient);
double gradient_th=sortedgradient[M-Mmax];
int m=0, i;
// Take the first Mused components with big gradients
for (i=0; i<M; i++)
if (optimizep[i] && Math.abs(gradient[i])>=gradient_th) {
m++; if (m==Mused) break;
} else optimizep[i]=false;
// Set the rest to 0
for (i=i+1; i<M; i++) optimizep[i]=false;
}
// Gradient descent
//for (int i=0; i<M; i++) if (optimizep[i]) x[i]-=0.01*gradient[i];
//if (true) return;
/* u will be a copy of the Hessian where we take only those
components corresponding to variables being optimized */
int kr=0, iw=0;
for (int ir = 0; ir<M; kr=kr+M,ir++) {
if (optimizep[ir]) {
int jw=0;
for (int jr = 0; jr<M; jr++)
if (optimizep[jr]) u[iw][jw++] = Hessian[kr + jr];
g[iw]=gradient[ir];
u[iw][iw] *= lambda;
iw++;
}
}
// Solve he equation system
/* SVD u=u*w*v^t */
update=unwarpJMathTools.linearLeastSquares(u,g);
/* x = x - update */
kr=0;
for (int kw = 0; kw<M; kw++)
if (optimizep[kw]) x[kw] -= update[kr++];
} /* end Marquardt_it */
/*--------------------------------------------------------------------------*/
private double optimizeCoeffs(
int intervals,
double thChangef,
double [][]cx,
double [][]cy
){
if (dialog!=null && dialog.isStopRegistrationSet()) return 0.0;
final double TINY = FLT_EPSILON;
final double EPS = 3.0e-8F;
final double FIRSTLAMBDA = 1;
final int MAXITER_OPTIMCOEFF = 300;
final int CUMULATIVE_SIZE = 5;
int int3=intervals+3;
int halfM=int3*int3;
int M=halfM*2;
double rescuedf, f;
double []x = new double [M];
double []rescuedx = new double [M];
double []diffx = new double [M];
double []rescuedgrad = new double [M];
double []grad = new double [M];
double []diffgrad = new double [M];
double []Hdx = new double [M];
double []rescuedhess = new double [M*M];
double []hess = new double [M*M];
double []safehess = new double [M*M];
double []proposedHess = new double [M*M];
boolean []optimize = new boolean [M];
int i, j, p, iter=1;
boolean skip_update, ill_hessian;
double improvementx=(double)Math.sqrt(TINY),
lambda=FIRSTLAMBDA, max_normx, distx, aux, gmax;
double fac, fae, dgdx, dxHdx, sumdiffg, sumdiffx;
unwarpJCumulativeQueue lastBest=
new unwarpJCumulativeQueue(CUMULATIVE_SIZE);
for (i=0; i<M; i++) optimize[i]=true;
/* Form the vector with the current guess for the optimization */
for (i= 0, p=0; i<intervals + 3; i++)
for (j = 0; j < intervals + 3; j++, p++) {
x[ p]=cx[i][j];
x[halfM+p]=cy[i][j];
}
/* Prepare the precomputed weights for interpolation */
swx = new unwarpJImageModel(x,intervals+3,intervals+3,0);
swy = new unwarpJImageModel(x,intervals+3,intervals+3,halfM);
swx.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
swy.precomputed_prepareForInterpolation(
target.getCurrentHeight(), target.getCurrentWidth(), intervals);
/* First computation of the similarity */
f=evaluateSimilarity(x,intervals,grad,false,false);
if (showMarquardtOptim) IJ.write("f(1)="+f);
/* Initially the hessian is the identity matrix multiplied by
the first function value */
for (i=0,p=0; i<M; i++)
for (j=0; j<M; j++,p++)
if (i==j) hess[p]=1.0F; else hess[p]=0.0F;
rescuedf = f;
for (i=0,p=0; i<M; i++) {
rescuedx[i]=x[i];
rescuedgrad[i]=grad[i];
for (j=0; j<M; j++,p++) rescuedhess[p]=hess[p];
}
int maxiter=MAXITER_OPTIMCOEFF*(source.getCurrentDepth()+1);
unwarpJProgressBar.stepProgressBar();
int last_successful_iter=0;
boolean stop=dialog!=null && dialog.isStopRegistrationSet();
while (iter < maxiter && !stop) {
/* Compute new x ------------------------------------------------- */
Marquardt_it(x, optimize, grad, hess, lambda);
/* Stopping criteria --------------------------------------------- */
/* Compute difference with the previous iteration */
max_normx=improvementx=0;
for (i=0; i<M; i++) {
diffx[i]=x[i]-rescuedx[i];
distx=Math.abs(diffx[i]);
improvementx+=distx*distx;
aux=Math.abs(rescuedx[i]) < Math.abs(x[i]) ? x[i] : rescuedx[i];
max_normx+=aux*aux;
}
if (TINY < max_normx) improvementx = improvementx/max_normx;
improvementx = (double) Math.sqrt(Math.sqrt(improvementx));
/* If there is no change with respect to the old geometry then
finish the iterations */
if (improvementx < Math.sqrt(TINY)) break;
/* Estimate the new function value -------------------------------- */
f=evaluateSimilarity(x,intervals,grad,false,false);
iter++;
if (showMarquardtOptim) IJ.write("f("+iter+")="+f+" lambda="+lambda);
unwarpJProgressBar.stepProgressBar();
/* Update lambda -------------------------------------------------- */
if (rescuedf > f) {
/* Check if the improvement is only residual */
lastBest.push_back(rescuedf-f);
if (lastBest.currentSize()==CUMULATIVE_SIZE && lastBest.getSum()/f<thChangef)
break;
/* If we have improved then estimate the hessian,
update the geometry, and decrease the lambda */
/* Estimate the hessian ....................................... */
if (showMarquardtOptim) IJ.write(" Accepted");
if ((last_successful_iter++%10)==0 && outputLevel>-1)
update_current_output(x,intervals);
/* Estimate the difference between gradients */
for (i=0; i<M; i++) diffgrad[i]=grad[i]-rescuedgrad[i];
/* Multiply this difference by the current inverse of the hessian */
for (i=0, p=0; i<M; i++) {
Hdx[i]=0.0F;
for (j=0; j<M; j++, p++) Hdx[i]+=hess[p]*diffx[j];
}
/* Calculate dot products for the denominators ................ */
dgdx=dxHdx=sumdiffg=sumdiffx=0.0F;
skip_update=true;
for (i=0; i<M; i++) {
dgdx += diffgrad[i]*diffx[i];
dxHdx += diffx[i]*Hdx[i];
sumdiffg += diffgrad[i]*diffgrad[i];
sumdiffx += diffx[i]*diffx[i];
if (Math.abs(grad[i])>=Math.abs(rescuedgrad[i])) gmax=Math.abs(grad[i]);
else gmax=Math.abs(rescuedgrad[i]);
if (gmax!=0 && Math.abs(diffgrad[i]-Hdx[i])>Math.sqrt(EPS)*gmax)
skip_update=false;
}
/* Update hessian ............................................. */
/* Skip if fac not sufficiently positive */
if (dgdx>Math.sqrt(EPS*sumdiffg*sumdiffx) && !skip_update) {
fae=1.0F/dxHdx;
fac=1.0F/dgdx;
/* Update the hessian after BFGS formula */
for (i=0, p=0; i<M; i++)
for (j=0; j<M; j++, p++) {
if (i<=j) proposedHess[p]=hess[p]+
fac*diffgrad[i]*diffgrad[j]
-fae*(Hdx[i]*Hdx[j]);
else proposedHess[p]=proposedHess[j*M+i];
}
ill_hessian=false;
if (!ill_hessian) {
for (i=0, p=0; i<M; i++)
for (j=0; j<M; j++,p++)
hess[p]= proposedHess[p];
} else
if (showMarquardtOptim)
IJ.write("Hessian cannot be safely updated, ill-conditioned");
} else
if (showMarquardtOptim)
IJ.write("Hessian cannot be safely updated");
/* Update geometry and lambda ................................. */
rescuedf = f;
for (i=0,p=0; i<M; i++) {
rescuedx[i]=x[i];
rescuedgrad[i]=grad[i];
for (j=0; j<M; j++,p++) rescuedhess[p]=hess[p];
}
if (1e-4 < lambda) lambda = lambda/10;
} else {
/* else, if it is worse, then recover the last geometry
and increase lambda, saturate lambda with FIRSTLAMBDA */
for (i=0,p=0; i<M; i++) {
x[i]=rescuedx[i];
grad[i]=rescuedgrad[i];
for (j=0; j<M; j++,p++) hess[p]=rescuedhess[p];
}
if (lambda < 1.0/TINY) lambda*=10;
else break;
if (lambda < FIRSTLAMBDA) lambda = FIRSTLAMBDA;
}
stop=dialog!=null && dialog.isStopRegistrationSet();
}
// Copy the values back to the input arrays
for (i= 0, p=0; i<intervals + 3; i++)
for (j = 0; j < intervals + 3; j++, p++) {
cx[i][j]=x[ p];
cy[i][j]=x[halfM+p];
}
unwarpJProgressBar.skipProgressBar(maxiter-iter);
return f;
}
/*-----------------------------------------------------------------------------*/
private double [][] propagateCoeffsToNextLevel(
int intervals,
final double [][]c,
double expansionFactor // Due to the change of size in the represented image
){
// Expand the coefficients for the next scale
intervals*=2;
double [][] cs_expand = new double[intervals+7][intervals+7];
// Upsample
for (int i=0; i<intervals+7; i++)
for (int j=0; j<intervals+7; j++) {
// If it is not in an even sample then set it to 0
if (i%2 ==0 || j%2 ==0) cs_expand[i][j]=0.0F;
else {
// Now look for this sample in the coarser level
int ipc=(i-1)/2;
int jpc=(j-1)/2;
cs_expand[i][j]=c[ipc][jpc];
}
}
// Define the FIR filter
double [][] u2n=new double [4][];
u2n[0]=null;
u2n[1]=new double[3]; u2n[1][0]=0.5F; u2n[1][1]=1.0F; u2n[1][2]=0.5F;
u2n[2]=null;
u2n[3]=new double[5]; u2n[3][0]=0.125F; u2n[3][1]=0.5F; u2n[3][2]=0.75F;
u2n[3][3]=0.5F; u2n[3][4]=0.125F;
int [] half_length_u2n={0, 1, 0, 2};
int kh=half_length_u2n[transformationSplineDegree];
// Apply the u2n filter to rows
double [][] cs_expand_aux = new double[intervals+7][intervals+7];
for (int i=1; i<intervals+7; i+=2)
for (int j=0; j<intervals+7; j++) {
cs_expand_aux[i][j]=0.0F;
for (int k=-kh; k<=kh; k++)
if (j+k>=0 && j+k<=intervals+6)
cs_expand_aux[i][j]+=u2n[transformationSplineDegree][k+kh]*cs_expand[i][j+k];
}
// Apply the u2n filter to columns
for (int i=0; i<intervals+7; i++)
for (int j=0; j<intervals+7; j++) {
cs_expand[i][j]=0.0F;
for (int k=-kh; k<=kh; k++)
if (i+k>=0 && i+k<=intervals+6)
cs_expand[i][j]+=u2n[transformationSplineDegree][k+kh]*cs_expand_aux[i+k][j];
}
// Copy the central coefficients to c
double [][]newc=new double [intervals+3][intervals+3];
for (int i= 0; i<intervals+3; i++)
for (int j = 0; j <intervals + 3; j++)
newc[i][j]=cs_expand[i+2][j+2]*expansionFactor;
// Return the new set of coefficients
return newc;
}
/*------------------------------------------------------------------*/
private void saveTransformation(
int intervals,
double [][]cx,
double [][]cy
) {
String filename=fn_tnf;
if (filename.equals("")) {
// Get the filename to save
File dir=new File(".");
String path="";
try {
path=dir.getCanonicalPath()+"/";
} catch (Exception e) {
e.printStackTrace();
}
filename=sourceImp.getTitle();
String new_filename="";
int dot = filename.lastIndexOf('.');
if (dot == -1) new_filename=filename + "_transf.txt";
else new_filename=filename.substring(0, dot)+"_transf.txt";
filename=path+filename;
if (outputLevel>-1) {
final Frame f = new Frame();
final FileDialog fd = new FileDialog(f, "Save Transformation", FileDialog.SAVE);
fd.setFile(new_filename);
fd.setVisible(true);
path = fd.getDirectory();
filename = fd.getFile();
if ((path == null) || (filename == null)) return;
filename=path+filename;
} else
filename=new_filename;
}
// Save the file
try {
final FileWriter fw = new FileWriter(filename);
String aux;
fw.write("Intervals="+intervals+"\n\n");
fw.write("X Coeffs -----------------------------------\n");
for (int i= 0; i<intervals + 3; i++) {
for (int j = 0; j < intervals + 3; j++) {
aux=""+cx[i][j];
while (aux.length()<21) aux=" "+aux;
fw.write(aux+" ");
}
fw.write("\n");
}
fw.write("\n");
fw.write("Y Coeffs -----------------------------------\n");
for (int i= 0; i<intervals + 3; i++) {
for (int j = 0; j < intervals + 3; j++) {
aux=""+cy[i][j];
while (aux.length()<21) aux=" "+aux;
fw.write(aux+" ");
}
fw.write("\n");
}
fw.close();
} catch (IOException e) {
IJ.error("IOException exception" + e);
} catch (SecurityException e) {
IJ.error("Security exception" + e);
}
}
/*------------------------------------------------------------------*/
private void showDeformationGrid(
int intervals,
double [][]cx,
double [][]cy,
ImageStack is2
) {
// Initialize output image
int stepv=Math.min(Math.max(10,targetCurrentHeight/15),30);
int stepu=Math.min(Math.max(10,targetCurrentWidth/15),30);
final double transformedImage [][]=new double [targetCurrentHeight][targetCurrentWidth];
for (int v=0; v<targetCurrentHeight; v++)
for (int u=0; u<targetCurrentWidth; u++) transformedImage[v][u]=255;
// Ask for memory for the transformation
double [][] transformation_x=new double [targetCurrentHeight][targetCurrentWidth];
double [][] transformation_y=new double [targetCurrentHeight][targetCurrentWidth];
// Compute the deformation
computeDeformation(intervals,cx,cy,transformation_x,transformation_y);
// Show deformed grid ........................................
// Show deformation vectors
for (int v=0; v<targetCurrentHeight; v+=stepv)
for (int u=0; u<targetCurrentWidth; u+=stepu) {
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
// Draw horizontal line
int uh=u+stepu;
if (uh<targetCurrentWidth) {
final double xh = transformation_x[v][uh];
final double yh = transformation_y[v][uh];
unwarpJMiscTools.drawLine(
transformedImage,
(int)Math.round(x) ,(int)Math.round(y),
(int)Math.round(xh),(int)Math.round(yh),0);
}
// Draw vertical line
int vv=v+stepv;
if (vv<targetCurrentHeight) {
final double xv = transformation_x[vv][u];
final double yv = transformation_y[vv][u];
unwarpJMiscTools.drawLine(
transformedImage,
(int)Math.round(x) ,(int)Math.round(y),
(int)Math.round(xv),(int)Math.round(yv),0);
}
}
// Set it to the image stack
FloatProcessor fp=new FloatProcessor(targetCurrentWidth,targetCurrentHeight);
for (int v=0; v<targetCurrentHeight; v++)
for (int u=0; u<targetCurrentWidth; u++)
fp.putPixelValue(u,v,transformedImage[v][u]);
is.addSlice("Deformation Grid",fp);
}
public ImageStack getIs() {
return is;
}
public void setIs(ImageStack is) {
this.is = is;
}
/*------------------------------------------------------------------*/
private void showDeformationVectors(
int intervals,
double [][]cx,
double [][]cy,
ImageStack is2
) {
// Initialize output image
int stepv=Math.min(Math.max(10,targetCurrentHeight/15),30);
int stepu=Math.min(Math.max(10,targetCurrentWidth/15),30);
final double transformedImage [][]=new double [targetCurrentHeight][targetCurrentWidth];
for (int v=0; v<targetCurrentHeight; v++)
for (int u=0; u<targetCurrentWidth; u++) transformedImage[v][u]=255;
// Ask for memory for the transformation
double [][] transformation_x=new double [targetCurrentHeight][targetCurrentWidth];
double [][] transformation_y=new double [targetCurrentHeight][targetCurrentWidth];
// Compute the deformation
computeDeformation(intervals,cx,cy,transformation_x,transformation_y);
// Show shift field ........................................
// Show deformation vectors
for (int v=0; v<targetCurrentHeight; v+=stepv)
for (int u=0; u<targetCurrentWidth; u+=stepu)
if (targetMsk.getValue(u,v)) {
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
if (sourceMsk.getValue(x,y))
unwarpJMiscTools.drawArrow(
transformedImage,
u,v,(int)Math.round(x),(int)Math.round(y),0,2);
}
// Set it to the image stack
FloatProcessor fp=new FloatProcessor(targetCurrentWidth,targetCurrentHeight);
for (int v=0; v<targetCurrentHeight; v++)
for (int u=0; u<targetCurrentWidth; u++)
fp.putPixelValue(u,v,transformedImage[v][u]);
is.addSlice("Deformation Field",fp);
}
/*-------------------------------------------------------------------*/
private void showTransformation(
final int intervals,
final double [][]cx, // Input, spline coefficients
final double [][]cy
){
boolean show_deformation=false;
// Ask for memory for the transformation
double [][] transformation_x=new double [targetHeight][targetWidth];
double [][] transformation_y=new double [targetHeight][targetWidth];
// Compute the deformation
computeDeformation(intervals,cx,cy,transformation_x,transformation_y);
if (show_deformation) {
unwarpJMiscTools.showImage("Transf. X",transformation_x);
unwarpJMiscTools.showImage("Transf. Y",transformation_y);
}
// Compute the warped image
FloatProcessor fp=new FloatProcessor(targetWidth,targetHeight);
FloatProcessor fp_mask=new FloatProcessor(targetWidth,targetHeight);
FloatProcessor fp_target=new FloatProcessor(targetWidth,targetHeight);
int uv=0;
for (int v=0; v<targetHeight; v++)
for (int u=0; u<targetWidth; u++,uv++) {
fp_target.putPixelValue(u,v,target.getImage()[uv]);
if (!targetMsk.getValue(u,v)) {
fp.putPixelValue(u,v,0);
fp_mask.putPixelValue(u,v,0);
} else {
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
if (sourceMsk.getValue(x,y)) {
source.prepareForInterpolation(x,y,ORIGINAL);
double sval=source.interpolateI();
fp.putPixelValue(u,v,sval);
fp_mask.putPixelValue(u,v,255);
} else {
fp.putPixelValue(u,v,0);
fp_mask.putPixelValue(u,v,0);
}
}
}
fp.resetMinAndMax();
//final ImageStack is = new ImageStack(targetWidth,targetHeight);
is = new ImageStack(targetWidth,targetHeight);
is.addSlice("Registered Image",fp);
System.out.println("outputLevel: " + outputLevel);
outputLevel = 2;
if (outputLevel>-1)
is.addSlice("Target Image",fp_target);
if (outputLevel>-1)
is.addSlice("Warped Source Mask",fp_mask);
if (outputLevel==2) {
System.out.println("true");
showDeformationVectors(intervals,cx,cy,is);
showDeformationGrid(intervals,cx,cy,is);
}
output_ip.setStack("Registered Image",is);
output_ip.setSlice(1); output_ip.getProcessor().resetMinAndMax();
output_ip.show();
if (outputLevel>-1) output_ip.updateAndRepaintWindow();
}
/*------------------------------------------------------------------*/
private void update_current_output(
final double []c,
int intervals
) {
// Set the coefficients to an interpolator
int cYdim=intervals+3;
int cXdim=cYdim;
int Nk=cYdim*cXdim;
swx.setCoefficients(c,cYdim,cXdim,0);
swy.setCoefficients(c,cYdim,cXdim,Nk);
// Compute the deformed image
FloatProcessor fp=new FloatProcessor(targetWidth,targetHeight);
int uv=0;
for (int v=0; v<targetHeight; v++)
for (int u=0; u<targetWidth; u++, uv++) {
if (targetMsk.getValue(u,v)) {
double down_u=u*factorWidth;
double down_v=v*factorHeight;
final double tv = (double)(down_v * intervals)/(double)(targetCurrentHeight-1) + 1.0F;
final double tu = (double)(down_u * intervals)/(double)(targetCurrentWidth -1) + 1.0F;
swx.prepareForInterpolation(tu,tv,ORIGINAL); double x=swx.interpolateI();
swy.prepareForInterpolation(tu,tv,ORIGINAL); double y=swy.interpolateI();
double up_x=x/factorWidth;
double up_y=y/factorHeight;
if (sourceMsk.getValue(up_x,up_y)) {
source.prepareForInterpolation(up_x,up_y,ORIGINAL);
fp.putPixelValue(u,v,target.getImage()[uv]-
source.interpolateI());
} else
fp.putPixelValue(u,v,0);
} else
fp.putPixelValue(u,v,0);
}
double min_val=output_ip.getProcessor().getMin();
double max_val=output_ip.getProcessor().getMax();
fp.setMinAndMax(min_val,max_val);
output_ip.setProcessor("Output",fp);
output_ip.updateImage();
// Draw the grid on the target image ...............................
// Some initialization
int stepv=Math.min(Math.max(10,targetHeight/15),30);
int stepu=Math.min(Math.max(10,targetWidth/15),30);
final double transformedImage [][]=new double [sourceHeight][sourceWidth];
double grid_colour=-1e-10;
uv=0;
for (int v=0; v<sourceHeight; v++)
for (int u=0; u<sourceWidth; u++,uv++) {
transformedImage[v][u]=source.getImage()[uv];
if (transformedImage[v][u]>grid_colour) grid_colour=transformedImage[v][u];
}
// Draw grid
for (int v=0; v<targetHeight+stepv; v+=stepv)
for (int u=0; u<targetWidth+stepu; u+=stepu) {
double down_u=u*factorWidth;
double down_v=v*factorHeight;
final double tv = (double)(down_v * intervals)/(double)(targetCurrentHeight-1) + 1.0F;
final double tu = (double)(down_u * intervals)/(double)(targetCurrentWidth -1) + 1.0F;
swx.prepareForInterpolation(tu,tv,ORIGINAL); double x=swx.interpolateI();
swy.prepareForInterpolation(tu,tv,ORIGINAL); double y=swy.interpolateI();
double up_x=x/factorWidth;
double up_y=y/factorHeight;
// Draw horizontal line
int uh=u+stepu;
if (uh<targetWidth+stepu) {
double down_uh=uh*factorWidth;
final double tuh = (double)(down_uh * intervals)/(double)(targetCurrentWidth -1) + 1.0F;
swx.prepareForInterpolation(tuh,tv,ORIGINAL); double xh=swx.interpolateI();
swy.prepareForInterpolation(tuh,tv,ORIGINAL); double yh=swy.interpolateI();
double up_xh=xh/factorWidth;
double up_yh=yh/factorHeight;
unwarpJMiscTools.drawLine(
transformedImage,
(int)Math.round(up_x) ,(int)Math.round(up_y),
(int)Math.round(up_xh),(int)Math.round(up_yh),grid_colour);
}
// Draw vertical line
int vv=v+stepv;
if (vv<targetHeight+stepv) {
double down_vv=vv*factorHeight;
final double tvv = (double)(down_vv * intervals)/(double)(targetCurrentHeight-1) + 1.0F;
swx.prepareForInterpolation(tu,tvv,ORIGINAL); double xv=swx.interpolateI();
swy.prepareForInterpolation(tu,tvv,ORIGINAL); double yv=swy.interpolateI();
double up_xv=xv/factorWidth;
double up_yv=yv/factorHeight;
unwarpJMiscTools.drawLine(
transformedImage,
(int)Math.round(up_x) ,(int)Math.round(up_y),
(int)Math.round(up_xv),(int)Math.round(up_yv),grid_colour);
}
}
// Update the target image plus
FloatProcessor fpg=new FloatProcessor(sourceWidth,sourceHeight);
for (int v=0; v<sourceHeight; v++)
for (int u=0; u<sourceWidth; u++)
fpg.putPixelValue(u,v,transformedImage[v][u]);
min_val=sourceImp.getProcessor().getMin();
max_val=sourceImp.getProcessor().getMax();
fpg.setMinAndMax(min_val,max_val);
sourceImp.setProcessor(sourceImp.getTitle(),fpg);
sourceImp.updateImage();
}
/*------------------------------------------------------------------*/
private double[] xWeight(
final double x,
final int xIntervals,
final boolean extended
) {
int length=xIntervals+1;
int j0=0, jF=xIntervals;
if (extended) {length+=2; j0--; jF++;}
final double[] b = new double[length];
final double interX = (double)xIntervals / (double)(targetCurrentWidth - 1);
for (int j=j0; j<=jF; j++) {
b[j-j0] = unwarpJMathTools.Bspline03(x * interX - (double)j);
}
return(b);
} /* end xWeight */
/*------------------------------------------------------------------*/
private double[] yWeight(
final double y,
final int yIntervals,
final boolean extended
) {
int length=yIntervals+1;
int i0=0, iF=yIntervals;
if (extended) {length+=2; i0--; iF++;}
final double[] b = new double[length];
final double interY = (double)yIntervals / (double)(targetCurrentHeight - 1);
for (int i = i0; i<=iF; i++) {
b[i-i0] = unwarpJMathTools.Bspline03(y * interY - (double)i);
}
return(b);
} /* end yWeight */
} /* end class unwarpJTransformation */