diff --git a/.gitignore b/.gitignore index e4e9e03d..11873ed5 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ log-*.txt.lck *.log *.gz package.sh -target/* \ No newline at end of file +target/* +/nbproject/ \ No newline at end of file diff --git a/src/main/java/ijfx/bridge/FxUserInterfaceBridge.java b/src/main/java/ijfx/bridge/FxUserInterfaceBridge.java index 17cf6ca4..765b6691 100644 --- a/src/main/java/ijfx/bridge/FxUserInterfaceBridge.java +++ b/src/main/java/ijfx/bridge/FxUserInterfaceBridge.java @@ -25,6 +25,7 @@ import ijfx.ui.datadisplay.image.ImageWindow; import ijfx.ui.main.ImageJFX; import ijfx.ui.datadisplay.table.TableWindow; +import ijfx.ui.datadisplay.text.TextWindow; import java.io.File; import java.util.concurrent.ExecutionException; import java.util.logging.Level; @@ -47,6 +48,7 @@ import org.scijava.ui.UserInterface; import org.scijava.ui.viewer.DisplayWindow; import mongis.utils.FXUtilities; +import org.scijava.display.TextDisplay; /** * UI bridge between ImageJFX and ImageJ @@ -204,7 +206,15 @@ public void show(final Display dspl) { ImageWindowContainer.getInstance().getChildren().add(new TableWindow((TableDisplay) dspl)); }); - } else { + } + else if (dspl instanceof TextDisplay) { + Platform.runLater(() -> { + + ImageWindowContainer.getInstance().getChildren().add(new TextWindow((TextDisplay) dspl)); + }); + + } + else { logger.warning("Cannot show display type :" + dspl.getClass().getSimpleName()); } } diff --git a/src/main/java/ijfx/plugins/AbstractImageJ1PluginAdapter.java b/src/main/java/ijfx/plugins/AbstractImageJ1PluginAdapter.java new file mode 100644 index 00000000..85c7d174 --- /dev/null +++ b/src/main/java/ijfx/plugins/AbstractImageJ1PluginAdapter.java @@ -0,0 +1,172 @@ +/* + * /* + * This file is part of ImageJ FX. + * + * ImageJ FX is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ImageJ FX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with ImageJ FX. If not, see . + * + * Copyright 2015,2016 Cyril MONGIS, Michael Knop + * + */ +package ijfx.plugins; + +import ij.ImagePlus; +import static java.lang.Math.toIntExact; +import java.util.stream.IntStream; +import net.imagej.Dataset; +import net.imagej.DatasetService; +import net.imagej.axis.Axes; +import net.imagej.axis.AxisType; +import net.imagej.axis.CalibratedAxis; +import net.imagej.display.ImageDisplay; +import net.imagej.display.ImageDisplayService; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.display.imagej.ImageJFunctions; + +import net.imglib2.type.numeric.integer.UnsignedShortType; +import org.scijava.ItemIO; +import org.scijava.command.Command; +import org.scijava.event.EventService; +import org.scijava.plugin.Parameter; + +/** + * + * @author Cyril MONGIS, 2015 + * @author Tuan anh TRINH + */ +public abstract class AbstractImageJ1PluginAdapter implements Command { + + @Parameter + DatasetService service; + + @Parameter + ImageDisplayService imageDisplayService; + + @Parameter + EventService eventService; + + @Parameter + boolean createCopy; + + boolean wholeWrap = false; + + public abstract ImagePlus processImagePlus(ImagePlus input); + + public ImagePlus getInput(Dataset dataset) { + return unwrapDataset(dataset); + } + + public Dataset setOutput(ImagePlus imp, Dataset dataset) { + if (!wholeWrap) { + dataset = wrapDataset(imp); + } else { + ImagePlus resultCopy = imp.duplicate(); + if (createCopy) { + dataset = emptyDataset(dataset, dataset.numDimensions()); + } + Dataset dataset2 = wrapDataset(resultCopy); + for (int i = 0; i < getNumberOfSlices(dataset); i++) { + dataset.setPlane(i, dataset2.getPlane(i)); + } + + } + return dataset; + } + + public static ImagePlus unwrapDataset(Dataset dataset) { + RandomAccessibleInterval r = (RandomAccessibleInterval) dataset.getImgPlus(); + ImagePlus wrapImage = ImageJFunctions.wrap(r, ""); + return wrapImage; + } + + public Dataset wrapDataset(ImagePlus imp) { + Img img = ImageJFunctions.wrap(imp.duplicate()); + return service.create(img); + } + + public static void configureImagePlus(ImagePlus imp, ImageDisplay imageDisplay) { + + imp.setC(imageDisplay.getIntPosition(Axes.CHANNEL)); + imp.setZ(imageDisplay.getIntPosition(Axes.Z)); + imp.setT(imageDisplay.getIntPosition(Axes.TIME)); + + } + + private Dataset emptyDataset(Dataset input, int sizeDims) { + AxisType[] axisType = new AxisType[input.numDimensions()]; + CalibratedAxis[] axeArray = new CalibratedAxis[input.numDimensions()]; + input.axes(axeArray); + + long[] dims = new long[sizeDims]; + for (int i = 0; i < sizeDims; i++) { + axisType[i] = axeArray[i].type(); + dims[i] = toIntExact(input.max(i) + 1); + } + return service.create(dims, input.getName(), axisType, input.getValidBits(), input.isSigned(), false); + } + + private Dataset chooseDataset(Dataset dataset) { + if (createCopy) { + return dataset.duplicateBlank(); + } + return dataset; + } + + public int getNumberOfSlices(Dataset dataset) { + + return (int) (dataset.getImgPlus().size() / (dataset.dimension(0) * dataset.dimension(1))); + + } + + public Dataset processDataset(Dataset dataset) { + Dataset datasetToModify = chooseDataset(dataset); + + for (int i = 0; i < getNumberOfSlices(dataset); i++) { + Dataset datasetOnePlane = emptyDataset(dataset, 2); + datasetOnePlane.setPlane(0, dataset.getPlane(i)); + ImagePlus result = processImagePlus(getInput(datasetOnePlane)); + setOutput(result.duplicate(), datasetOnePlane); + datasetToModify.setPlane(i, datasetOnePlane.getPlane(0)); + } + dataset = datasetToModify; + return dataset; + } +// +// /** +// * Wrap the whole Dataset. Use more memory +// * +// * @param dataset +// * @return +// */ +// public Dataset processDatasetWholeWrap(Dataset dataset) { +// ImagePlus result = processImagePlus(getInput(dataset)); +// ImagePlus resultCopy = result.duplicate(); +// if (createCopy) { +// dataset = emptyDataset(dataset, dataset.numDimensions()); +// } +// Dataset dataset2 = wrapDataset(resultCopy); +// for (int i = 0; i < getNumberOfSlices(dataset); i++) { +// dataset.setPlane(i, dataset2.getPlane(i)); +// } +// return dataset; +// } + + public boolean isWholeWrap() { + return wholeWrap; + } + + public void setWholeWrap(boolean wholeWrap) { + this.wholeWrap = wholeWrap; + } +} diff --git a/src/main/java/ijfx/plugins/DefaultImageJ1PluginAdapter.java b/src/main/java/ijfx/plugins/DefaultImageJ1PluginAdapter.java new file mode 100644 index 00000000..fe8b66fc --- /dev/null +++ b/src/main/java/ijfx/plugins/DefaultImageJ1PluginAdapter.java @@ -0,0 +1,49 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins; + +import ij.ImagePlus; +import net.imagej.Dataset; +import org.scijava.ItemIO; +import org.scijava.command.Command; +import org.scijava.plugin.Attr; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +/** + * + * @author Tuan anh TRINH + */ +@Plugin(type = Command.class, menuPath = "Plugins>DefaultAdapter") +public class DefaultImageJ1PluginAdapter extends AbstractImageJ1PluginAdapter { + @Parameter(type = ItemIO.BOTH) + protected Dataset dataset; + @Override + public ImagePlus processImagePlus(ImagePlus input) { + return input; + } + + @Override + public void run() { + dataset = processDataset(dataset); + } + + + +} diff --git a/src/main/java/ijfx/plugins/DefaultWholeWrapper.java b/src/main/java/ijfx/plugins/DefaultWholeWrapper.java new file mode 100644 index 00000000..b35c60f3 --- /dev/null +++ b/src/main/java/ijfx/plugins/DefaultWholeWrapper.java @@ -0,0 +1,51 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins; + +import ij.ImagePlus; +import java.io.File; +import net.imagej.Dataset; +import org.scijava.ItemIO; +import org.scijava.command.Command; +import org.scijava.plugin.Attr; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +/** + * + * @author Tuan anh TRINH + */ +@Plugin(type = Command.class, menuPath = "Plugins>DefaultWholeAdapter") +public class DefaultWholeWrapper extends AbstractImageJ1PluginAdapter { + @Parameter(type = ItemIO.BOTH) + protected Dataset dataset; + @Parameter + File f; + @Override + public ImagePlus processImagePlus(ImagePlus input) { + return input; + } + + @Override + public void run() { + setWholeWrap(true); + dataset = setOutput(getInput(dataset), dataset); + } + +} diff --git a/src/main/java/ijfx/plugins/UnwarpJ/UnwarpJ_.java b/src/main/java/ijfx/plugins/UnwarpJ/UnwarpJ_.java new file mode 100644 index 00000000..91426a2b --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/UnwarpJ_.java @@ -0,0 +1,568 @@ +/* + * /* + * This file is part of ImageJ FX. + * + * ImageJ FX is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ImageJ FX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with ImageJ FX. If not, see . + * + * Copyright 2015,2016 Cyril MONGIS, Michael Knop + * + */ +/*==================================================================== +| Version: July 8, 2006 by Steve Bryson and Carlos Oscar +| 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. +\===================================================================*/ + +package ijfx.plugins.UnwarpJ; + +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.plugin.JpegWriter; +import ij.io.Opener; +import ij.measure.Calibration; +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 ijfx.plugins.AbstractImageJ1PluginAdapter; +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; +import net.imagej.Dataset; +import org.scijava.ItemIO; +import org.scijava.command.Command; +import org.scijava.plugin.Attr; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; + +/*==================================================================== +| UnwarpJ_ +\===================================================================*/ +/* Main class. + This is the one called by ImageJ */ + +/*------------------------------------------------------------------*/ + +@Plugin(type = Command.class, menuPath = "Plugins>UnwarpJ") +public class UnwarpJ_ extends AbstractImageJ1PluginAdapter { + @Parameter + protected Dataset dataset; + @Parameter(label = "Source") + Dataset source; + + @Parameter(label = "Target") + Dataset target; + /* begin class UnwarpJ_ */ + +/*.................................................................... + Public methods +....................................................................*/ + +/*------------------------------------------------------------------*/ + @Override + public void run() + { + setWholeWrap(true); + 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[]) { +String []args= new String("unwarpj -align target.jpg NULL source.jpg NULL 0 2 1 1 0 output.tif 1 landmarks.txt").split(" "); + + 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 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; + 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 = this.getInput(this.dataset); +// 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 = getInput(dataset); +// 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); +} + + @Override + public ImagePlus processImagePlus(ImagePlus input) { + return input; + } + + + +} /* end class UnwarpJ_ */ + +/*==================================================================== +| unwarpJClearAll +\===================================================================*/ + + +/*==================================================================== +| unwarpJCredits +\===================================================================*/ + + + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ + + +/*==================================================================== +| unwarpJFile +\===================================================================*/ + + + +/*==================================================================== +| unwarpJFinalAction +\===================================================================*/ + + +/*==================================================================== +| unwarpJImageModel +\===================================================================*/ + + +/*==================================================================== +| unwarpJMask +\===================================================================*/ + + + + +/*==================================================================== +| unwarpJPointAction +\===================================================================*/ + + +/*==================================================================== +| unwarpJPointHandler +\===================================================================*/ + + +/*==================================================================== +| unwarpJPointToolbar +\===================================================================*/ + + +/*==================================================================== +| unwarpJProgressBar +\===================================================================*/ + + +/*==================================================================== +| unwarpJTransformation +\===================================================================*/ + diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJClearAll.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJClearAll.java new file mode 100644 index 00000000..e95e62de --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJClearAll.java @@ -0,0 +1,87 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.ImagePlus; +import java.awt.Button; +import java.awt.Dialog; +import java.awt.Frame; +import java.awt.GridLayout; +import java.awt.Insets; +import java.awt.Label; +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; + +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCredits.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCredits.java new file mode 100644 index 00000000..e22d4098 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCredits.java @@ -0,0 +1,94 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import java.awt.BorderLayout; +import java.awt.Button; +import java.awt.Dialog; +import java.awt.FlowLayout; +import java.awt.Frame; +import java.awt.Insets; +import java.awt.Label; +import java.awt.Panel; +import java.awt.TextArea; +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; + +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCumulativeQueue.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCumulativeQueue.java new file mode 100644 index 00000000..6a963372 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJCumulativeQueue.java @@ -0,0 +1,82 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import java.util.Vector; + +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJDialog.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJDialog.java new file mode 100644 index 00000000..b234bf08 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJDialog.java @@ -0,0 +1,796 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.gui.GUI; +import ij.gui.ImageCanvas; +import ij.gui.Toolbar; +import ij.process.FloatProcessor; +import java.awt.Button; +import java.awt.Checkbox; +import java.awt.Choice; +import java.awt.Dialog; +import java.awt.FlowLayout; +import java.awt.Frame; +import java.awt.GridLayout; +import java.awt.Insets; +import java.awt.Label; +import java.awt.Panel; +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.TextEvent; +import java.awt.event.TextListener; + +/*==================================================================== +| 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(); + 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFile.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFile.java new file mode 100644 index 00000000..9766e365 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFile.java @@ -0,0 +1,263 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import java.awt.Button; +import java.awt.CheckboxGroup; +import java.awt.Dialog; +import java.awt.FileDialog; +import java.awt.Font; +import java.awt.Frame; +import java.awt.GridLayout; +import java.awt.Insets; +import java.awt.Label; +import java.awt.Point; +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; +import java.io.FileWriter; +import java.io.IOException; +import java.util.Stack; +import java.util.Vector; + +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFinalAction.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFinalAction.java new file mode 100644 index 00000000..3edad90b --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJFinalAction.java @@ -0,0 +1,147 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.ImagePlus; +import ij.gui.ImageWindow; +import ij.process.FloatProcessor; + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ +/*==================================================================== +| 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 unwarpJFinalAction + * 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 unwarpJDialog to + * unwarpJFinalAction. + ********************************************************************/ + 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 unwarpJDialog 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); + } + +} diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJImageModel.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJImageModel.java new file mode 100644 index 00000000..10c8d629 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJImageModel.java @@ -0,0 +1,1329 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.process.ImageProcessor; +import java.util.Stack; + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ +/*==================================================================== +| unwarpJFile +\===================================================================*/ +/*==================================================================== +| unwarpJFinalAction +\===================================================================*/ +/*==================================================================== +| 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 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; + double 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; + double sxx = 0.0F; + double 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 + ) { + // 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 + ) { + // 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMask.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMask.java new file mode 100644 index 00000000..850054c7 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMask.java @@ -0,0 +1,238 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.process.ByteProcessor; +import ij.process.FloatProcessor; +import ij.process.ImageProcessor; +import ij.process.ShortProcessor; +import java.awt.Point; +import java.awt.Polygon; +import java.awt.Rectangle; +import java.util.Vector; + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ +/*==================================================================== +| unwarpJFile +\===================================================================*/ +/*==================================================================== +| unwarpJFinalAction +\===================================================================*/ +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMathTools.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMathTools.java new file mode 100644 index 00000000..18f69581 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMathTools.java @@ -0,0 +1,547 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ +/*==================================================================== +| unwarpJFile +\===================================================================*/ +/*==================================================================== +| unwarpJFinalAction +\===================================================================*/ +/*==================================================================== +| unwarpJImageModel +\===================================================================*/ +/*==================================================================== +| 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; + double scale; + double c; + double f; + double g; + double h; + double s; + double x; + double y; + double 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMiscTools.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMiscTools.java new file mode 100644 index 00000000..caf00174 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJMiscTools.java @@ -0,0 +1,574 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.gui.ImageWindow; +import ij.process.ByteProcessor; +import ij.process.FloatProcessor; +import ij.process.ImageProcessor; +import ij.process.ShortProcessor; +import java.io.BufferedReader; +import java.io.FileNotFoundException; +import java.io.FileReader; +import java.io.IOException; +import java.util.Stack; +import java.util.StringTokenizer; + +/*==================================================================== +| unwarpJDialog +\===================================================================*/ +/*==================================================================== +| unwarpJFile +\===================================================================*/ +/*==================================================================== +| unwarpJFinalAction +\===================================================================*/ +/*==================================================================== +| unwarpJImageModel +\===================================================================*/ +/*==================================================================== +| unwarpJMask +\===================================================================*/ +/*==================================================================== +| 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. */ + public static 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) */ + public static 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); + } + } 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); + } + } 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 + java.awt.Point p2 = new java.awt.Point(-arrow_size2, -arrow_size); + java.awt.Point p3 = new java.awt.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 */ + public static 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; + int 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; + } + } + + /*------------------------------------------------------------------*/ + public static 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]; + } + } + } + + public static 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. */ + public static void loadPoints(String filename, Stack sourceStack, Stack targetStack) { + java.awt.Point sourcePoint; + java.awt.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# #xSource# #ySource# #xTarget# #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# #xSource# #ySource# #xTarget# #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# #xSource# #ySource# #xTarget# #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# #xSource# #ySource# #xTarget# #yTarget#"); + return; + } + xTarget = line.substring(0, separatorIndex); + xTarget = xTarget.trim(); + yTarget = line.substring(separatorIndex); + yTarget = yTarget.trim(); + sourcePoint = new java.awt.Point(Integer.valueOf(xSource).intValue(), Integer.valueOf(ySource).intValue()); + sourceStack.push(sourcePoint); + targetPoint = new java.awt.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; + } + } + + public static 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; + } + } + + /*------------------------------------------------------------------*/ + public static 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. */ + public static 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointAction.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointAction.java new file mode 100644 index 00000000..992995af --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointAction.java @@ -0,0 +1,295 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.WindowManager; +import ij.gui.ImageCanvas; +import ij.measure.Calibration; +import java.awt.Event; +import java.awt.Point; +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; + +/*==================================================================== +| 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(); + int xp; + int y = e.getY(); + int 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointHandler.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointHandler.java new file mode 100644 index 00000000..b28790f7 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointHandler.java @@ -0,0 +1,455 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.WindowManager; +import ij.gui.ImageCanvas; +import ij.gui.ImageWindow; +import ij.gui.Roi; +import java.awt.Color; +import java.awt.Graphics; +import java.awt.Point; +import java.util.Vector; + +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointToolbar.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointToolbar.java new file mode 100644 index 00000000..519a3df5 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJPointToolbar.java @@ -0,0 +1,524 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.gui.GUI; +import ij.gui.Toolbar; +import java.awt.Canvas; +import java.awt.Color; +import java.awt.Component; +import java.awt.Container; +import java.awt.Graphics; +import java.awt.event.MouseEvent; +import java.awt.event.MouseListener; + +/*==================================================================== +| unwarpJPointAction +\===================================================================*/ +/*==================================================================== +| 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJProgressBar.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJProgressBar.java new file mode 100644 index 00000000..1a1b9623 --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJProgressBar.java @@ -0,0 +1,116 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; + +/*==================================================================== +| unwarpJPointAction +\===================================================================*/ +/*==================================================================== +| unwarpJPointHandler +\===================================================================*/ +/*==================================================================== +| 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 batch. + * + * @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 stride 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 1 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 batch 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJ/unwarpJTransformation.java b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJTransformation.java new file mode 100644 index 00000000..5b58a1bf --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJ/unwarpJTransformation.java @@ -0,0 +1,2260 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins.UnwarpJ; + +import ij.IJ; +import ij.ImagePlus; +import ij.ImageStack; +import ij.process.FloatProcessor; +import java.awt.FileDialog; +import java.awt.Frame; +import java.awt.Point; +import java.io.File; +import java.io.FileWriter; +import java.io.IOException; +import java.util.Arrays; +import java.util.Vector; + +/*==================================================================== +| 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 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 + ) { + 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 + ) { + 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 + ) { + // 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]; + double 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; + double gradcurly = 0.0; + double graddivx = 0.0; + double graddivy = 0.0; + double xdx = 0.0; + double xdy = 0.0; + double ydx = 0.0; + double ydy = 0.0; + double xdxdy = 0.0; + double xdxdx = 0.0; + double xdydy = 0.0; + double ydxdy = 0.0; + double ydxdx = 0.0; + double 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; + int 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 -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 is) { + // 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); + } + + /*------------------------------------------------------------------*/ + private void showDeformationVectors(int intervals, double[][] cx, double[][] cy, ImageStack is) { + // 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.addSlice("Registered Image", fp); + if (outputLevel > -1) { + is.addSlice("Target Image", fp_target); + } + if (outputLevel > -1) { + is.addSlice("Warped Source Mask", fp_mask); + } + if (outputLevel == 2) { + showDeformationVectors(intervals, cx, cy, is); + showDeformationGrid(intervals, cx, cy, is); + } + output_ip.setStack("Registered Image", is); + output_ip.setSlice(1); + output_ip.getProcessor().resetMinAndMax(); + 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; + int 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; + int 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 */ \ No newline at end of file diff --git a/src/main/java/ijfx/plugins/UnwarpJPlugin.java b/src/main/java/ijfx/plugins/UnwarpJPlugin.java new file mode 100644 index 00000000..ea8cae2f --- /dev/null +++ b/src/main/java/ijfx/plugins/UnwarpJPlugin.java @@ -0,0 +1,49 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.plugins; + +import ij.ImagePlus; +import net.imagej.Dataset; +import org.scijava.command.Command; +import org.scijava.plugin.Attr; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; + +@Plugin(type = Command.class, menuPath = "Plugins>UnwarpJ", attrs = { + @Attr(name = "no-legacy")}) +public class UnwarpJPlugin extends AbstractImageJ1PluginAdapter { + + @Parameter(label = "Source") + Dataset source; + + @Parameter(label = "Target") + Dataset target; + + @Override + public ImagePlus processImagePlus(ImagePlus input) { + return input; + } + + @Override + public void run() { + throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates. + } + +} diff --git a/src/main/java/ijfx/plugins/commands/MagicWand.java b/src/main/java/ijfx/plugins/commands/MagicWand.java index 913c4caa..47260bf7 100644 --- a/src/main/java/ijfx/plugins/commands/MagicWand.java +++ b/src/main/java/ijfx/plugins/commands/MagicWand.java @@ -26,7 +26,7 @@ import ij.plugin.frame.Recorder; import ij.process.ByteProcessor; import ij.process.ImageProcessor; -import ijfx.plugins.ImageJ1PluginAdapter; +import ijfx.plugins.AbstractImageJ1PluginAdapter; import ijfx.service.object_detection.ObjectDetectionService; import java.util.List; import net.imagej.Dataset; @@ -98,9 +98,9 @@ public void run() { getPixelsNears(pos, dataset.randomAccess()); - ImagePlus imp = ImageJ1PluginAdapter.unwrapDataset(dataset); + ImagePlus imp = AbstractImageJ1PluginAdapter.unwrapDataset(dataset); - ImageJ1PluginAdapter.configureImagePlus(imp, imageDisplay); + AbstractImageJ1PluginAdapter.configureImagePlus(imp, imageDisplay); doWand(imp, originX, originY, 1); diff --git a/src/main/java/ijfx/service/preview/PreviewService.java b/src/main/java/ijfx/service/preview/PreviewService.java index 0033fe02..b1f57498 100644 --- a/src/main/java/ijfx/service/preview/PreviewService.java +++ b/src/main/java/ijfx/service/preview/PreviewService.java @@ -23,6 +23,7 @@ import ijfx.service.batch.BatchSingleInput; import ijfx.service.batch.DisplayBatchInput; import java.awt.image.BufferedImage; +import static java.lang.Math.abs; import java.util.List; import java.util.Map; import java.util.logging.Logger; @@ -74,7 +75,7 @@ public class PreviewService extends AbstractService implements ImageJService { @Parameter ModuleService moduleService; - + @Parameter BatchService batchService; private int width; @@ -119,11 +120,17 @@ public void setParameters(int x, int y, int w, int h) { this.y = y; this.width = w; this.height = h; + int widthDataset = (int) imageDisplayService.getActiveDataset().max(0); + int heightDataset = (int) imageDisplayService.getActiveDataset().max(1); if (x < 0 || y < 0) { - int widthDataset = (int) imageDisplayService.getActiveDataset().max(0); - int heightDataset = (int) imageDisplayService.getActiveDataset().max(1); this.x = (int) (widthDataset / 2.0 - width / 2.0); this.y = (int) (heightDataset / 2.0 - height / 2.0); + } + if (widthDataset < width || heightDataset < height) { + this.x = 0; + this.y = 0; + width = widthDataset; + height = heightDataset; } } @@ -149,15 +156,13 @@ public Image getImageDisplay(String command, Map inputMap) { * @return output */ public Dataset getEmptyDataset(Dataset input) { - AxisType[] axisType = new AxisType[input.numDimensions()]; - CalibratedAxis[] axeArray = new CalibratedAxis[input.numDimensions()]; + AxisType[] axisType = new AxisType[2]; + CalibratedAxis[] axeArray = new CalibratedAxis[2]; input.axes(axeArray); - long[] dims = new long[axeArray.length]; + long[] dims = new long[2]; for (int i = 0; i < dims.length; i++) { axisType[i] = axeArray[i].type(); - dims[i] = 1;// toIntExact(input.max(i) + 1); - } dims[0] = width; dims[1] = height; @@ -259,11 +264,16 @@ public Dataset applyCommand(Dataset dataset, String command, Map BatchSingleInput batchSingleInput = new DisplayBatchInput(); this.context().inject(batchSingleInput); batchSingleInput.setDataset(dataset); - + Module module = moduleService.createModule(commandService.getCommand(command)); //CommandInfo commandInfo = new CommandInfo(command); - // Module module = new CommandModule(commandInfo); - this.context().inject(module.getDelegateObject()); + // Module module = new CommandModule(commandInfo); + try { + this.context().inject(module.getDelegateObject()); + + } catch (Exception e) { + e.printStackTrace(); + } batchService.executeModule(batchSingleInput, module, inputMap); Dataset result = batchSingleInput.getDataset(); return result; diff --git a/src/main/java/ijfx/ui/datadisplay/text/TextDisplayView.java b/src/main/java/ijfx/ui/datadisplay/text/TextDisplayView.java index 914c7c82..fa821ac3 100644 --- a/src/main/java/ijfx/ui/datadisplay/text/TextDisplayView.java +++ b/src/main/java/ijfx/ui/datadisplay/text/TextDisplayView.java @@ -20,10 +20,50 @@ */ package ijfx.ui.datadisplay.text; +import ijfx.ui.main.ImageJFX; +import java.io.IOException; +import java.util.logging.Level; +import java.util.logging.Logger; +import javafx.fxml.FXML; +import javafx.scene.layout.BorderPane; +import javafx.scene.text.Text; +import mongis.utils.FXUtilities; +import org.scijava.display.TextDisplay; + /** * * @author Cyril MONGIS, 2015 */ -public class TextDisplayView { +public class TextDisplayView extends BorderPane { + + @FXML + Text text; + + TextDisplay textDisplay; + + + final Logger logger = ImageJFX.getLogger(); + + + public TextDisplayView() { + + logger.info("Injecting FXML"); + try { + // inject TextDisplayView.fxml from the class name + FXUtilities.injectFXML(this, "/ijfx/ui/text/TextDisplayView.fxml"); + logger.info("FXML injected"); + } catch (IOException ex) { + ImageJFX.getLogger().log(Level.SEVERE, null, ex); + } + logger.info("Creating text model"); + logger.info("Text created"); + + } + + TextDisplayView(TextDisplay textDisplay) { + this(); + this.textDisplay = textDisplay; + text.setText(this.textDisplay.get(0)); + } } diff --git a/src/main/java/ijfx/ui/datadisplay/text/TextWindow.java b/src/main/java/ijfx/ui/datadisplay/text/TextWindow.java new file mode 100644 index 00000000..91d960f6 --- /dev/null +++ b/src/main/java/ijfx/ui/datadisplay/text/TextWindow.java @@ -0,0 +1,96 @@ +/* + This file is part of ImageJ FX. + + ImageJ FX is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ImageJ FX is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ImageJ FX. If not, see . + + Copyright 2015,2016 Cyril MONGIS, Michael Knop + + */ +package ijfx.ui.datadisplay.text; + +import ijfx.ui.main.ImageJFX; +import java.util.logging.Logger; +import javafx.event.EventType; +import javafx.scene.control.Button; +import javafx.scene.input.MouseEvent; +import javafx.scene.layout.BorderPane; +import jfxtras.scene.control.window.CloseIcon; +import jfxtras.scene.control.window.Window; +import org.scijava.display.DisplayService; +import org.scijava.display.TextDisplay; +import org.scijava.display.event.DisplayActivatedEvent; +import org.scijava.display.event.DisplayCreatedEvent; +import org.scijava.display.event.DisplayDeletedEvent; +import org.scijava.event.EventHandler; +import org.scijava.event.EventService; +import org.scijava.plugin.Parameter; + +/** + * + * @author Tuan anh TRINH + */ +public class TextWindow extends Window { + + Logger logger = ImageJFX.getLogger(); + @Parameter + DisplayService displayService; + + @Parameter + EventService eventService; + private TextDisplay textDisplay; + + public TextWindow() { + super(); + CloseIcon closeIcon = new CloseIcon(this); + getRightIcons().add(closeIcon); + + } + + public TextWindow(TextDisplay display) { + this(); + this.setTitle(display.getName()); + textDisplay = display; + setContentPane(new TextDisplayView(textDisplay)); + for (EventType t : new EventType[]{MouseEvent.MOUSE_CLICKED, MouseEvent.DRAG_DETECTED, MouseEvent.MOUSE_PRESSED}) { + addEventHandler(t, this::putInFront); + getContentPane().addEventHandler(t, this::putInFront); + } + textDisplay.getContext().inject(this); + this.setOnCloseAction(event -> { + eventService.publishLater(new DisplayDeletedEvent(textDisplay)); + }); + + } + + public void putInFront(MouseEvent event) { + putInFront(); + } + + public void putInFront() { + logger.info("Putting in front " + textDisplay); + displayService.setActiveDisplay(textDisplay); + displayService.getContext().toString(); + setMoveToFront(true); + } + + @EventHandler + protected void onActiveDisplayChanged(DisplayActivatedEvent event) { + if (event.getDisplay() == textDisplay) { + System.out.println("I'm activated !" + textDisplay.getName()); + setFocused(true); + } else { + setFocused(false); + } + } +} diff --git a/src/main/java/ijfx/ui/explorer/DefaultFolder.java b/src/main/java/ijfx/ui/explorer/DefaultFolder.java index c8e6485f..4edf2de1 100644 --- a/src/main/java/ijfx/ui/explorer/DefaultFolder.java +++ b/src/main/java/ijfx/ui/explorer/DefaultFolder.java @@ -222,28 +222,7 @@ private void addItems(List explorables) { } - /* - private Integer fetchMoreStatistics(List explorableList) { - Integer elementAnalyzedCount = 0; - int elements = explorableList.size(); - int i = 0; - for (Explorable e : explorableList) { - statusService.showStatus(i, elements, "Fetchting min/max for more exploration."); - - if (!e.getMetaDataSet().containsKey(MetaData.STATS_PIXEL_MIN)) { - if (e instanceof ImageRecordIconizer) { - ImageRecordIconizer iconizer = (ImageRecordIconizer) e; - SummaryStatistics stats = statsService.getStatistics(iconizer.getImageRecord().getFile()); - iconizer.getMetaDataSet().putGeneric(MetaData.STATS_PIXEL_MIN, stats.getMin()); - iconizer.getMetaDataSet().putGeneric(MetaData.STATS_PIXEL_MAX, stats.getMax()); - iconizer.getMetaDataSet().putGeneric(MetaData.STATS_PIXEL_MEAN, stats.getMean()); - elementAnalyzedCount++; - } - } - i++; - } - return elementAnalyzedCount; - }*/ + diff --git a/src/main/java/ijfx/ui/explorer/view/GridIconView.java b/src/main/java/ijfx/ui/explorer/view/GridIconView.java index ee3fea30..28de1e76 100644 --- a/src/main/java/ijfx/ui/explorer/view/GridIconView.java +++ b/src/main/java/ijfx/ui/explorer/view/GridIconView.java @@ -73,8 +73,8 @@ public GridIconView() { topBar = new GridPane(); listLabel = new ArrayList<>(); - listLabel.add(new Label("Columns")); listLabel.add(new Label("Rows")); + listLabel.add(new Label("Columns")); listLabel.add(new Label("Group by")); IntStream.range(0, listLabel.size()).forEach(i -> { topBar.add(listLabel.get(i), i, 0); diff --git a/src/main/java/ijfx/ui/explorer/view/GroupExplorable.java b/src/main/java/ijfx/ui/explorer/view/GroupExplorable.java index 3109026e..eb310418 100644 --- a/src/main/java/ijfx/ui/explorer/view/GroupExplorable.java +++ b/src/main/java/ijfx/ui/explorer/view/GroupExplorable.java @@ -85,7 +85,7 @@ public boolean checkNumber(String metaData){ { List filtered = filterExplorable(listItems, metaData); SortExplorableUtils.sort(metaData, filtered); - return SortExplorableUtils.findLimits(metaData, filtered).size() <= 10; + return SortExplorableUtils.findLimits(metaData, filtered).size() <= 25; } } diff --git a/src/main/java/ijfx/ui/module/skin/DatasetComboBoxInputSkin.java b/src/main/java/ijfx/ui/module/skin/DatasetComboBoxInputSkin.java index 98ec2bf3..be1f9a47 100644 --- a/src/main/java/ijfx/ui/module/skin/DatasetComboBoxInputSkin.java +++ b/src/main/java/ijfx/ui/module/skin/DatasetComboBoxInputSkin.java @@ -24,7 +24,9 @@ import ijfx.ui.module.input.Input; import java.util.List; import java.util.stream.Collectors; +import javafx.beans.property.ObjectProperty; import javafx.beans.property.Property; +import javafx.beans.property.SimpleObjectProperty; import javafx.scene.Node; import javafx.scene.control.ComboBox; import net.imagej.Dataset; @@ -37,11 +39,12 @@ /** * * @author Cyril MONGIS, 2015 + * @author Tuan anh TRINH */ @Plugin(type = InputSkinPlugin.class) public class DatasetComboBoxInputSkin extends AbstractInputSkinPlugin { - //ObjectProperty valueProperty = new SimpleObjectProperty<>(); + ObjectProperty valueProperty = new SimpleObjectProperty<>(); @Parameter ImageDisplayService imageDisplayService; @@ -56,7 +59,7 @@ public class DatasetComboBoxInputSkin extends AbstractInputSkinPlugin { @Override public Property valueProperty() { - return datasetComboBox.itemsProperty(); + return valueProperty; } @Override @@ -71,25 +74,24 @@ public void dispose() { @Override public boolean canHandle(Class clazz) { - //return false; - return clazz.isAssignableFrom(Dataset.class); + System.out.println("can i handle ?" + clazz); + return clazz == Dataset.class; } @Override public void init(Input input) { //datasetComboBox.getItems().addAll(datasetService.getDatasets()); - - - List toAdd = datasetService.getDatasets() + List toAdd = datasetService.getDatasets() .parallelStream() .filter(dataset -> dataset.toString().endsWith("lut") == false) .collect(Collectors.toList()); - datasetComboBox.getItems().addAll(toAdd); - + datasetComboBox.getItems().addAll(toAdd); + datasetComboBox.getSelectionModel().selectedItemProperty().addListener((obs, oldValue, newValue) -> { + valueProperty.setValue(newValue); + }); - //input.getValue(); } diff --git a/src/main/resources/ijfx/ui/main/flatterfx.css b/src/main/resources/ijfx/ui/main/flatterfx.css index 53e3cd5b..5035fbbf 100644 --- a/src/main/resources/ijfx/ui/main/flatterfx.css +++ b/src/main/resources/ijfx/ui/main/flatterfx.css @@ -1928,4 +1928,4 @@ list-cell:filled:selected:focused, .list-cell:filled:selected{ -glyph-size:150; -fx-translate-x:10; -glyph-color:red; - } \ No newline at end of file + } diff --git a/src/main/resources/ijfx/ui/menutoolbar/toolbarSettings.json b/src/main/resources/ijfx/ui/menutoolbar/toolbarSettings.json index 468fded9..1b35eb87 100644 --- a/src/main/resources/ijfx/ui/menutoolbar/toolbarSettings.json +++ b/src/main/resources/ijfx/ui/menutoolbar/toolbarSettings.json @@ -513,4 +513,4 @@ -] \ No newline at end of file +] diff --git a/src/main/resources/ijfx/ui/text/TextDisplayView.fxml b/src/main/resources/ijfx/ui/text/TextDisplayView.fxml new file mode 100644 index 00000000..ceacccc1 --- /dev/null +++ b/src/main/resources/ijfx/ui/text/TextDisplayView.fxml @@ -0,0 +1,20 @@ + + + + + + + + + + + + +
+ + + + + +
+