Java tutorial
/** ImageJ plugin Copyright (C) 2008 Verena Kaynig. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt ) This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. **/ /* **************************************************************************** * * This ImageJ Plugin estimates a non linear lens distortion in * * the images given and corrects all images with the same transformation * * * * * * TODO: * * - show histogram of beta coefficients to evaluate kernel dimension needed * * - show SIFT matches before and after correction? * * - give numerical xcorr results in a better manner * * - show visual impression of the stitching ? * * - DOKUMENTATION * * - do mask images properly to incorporate into TrakEM * * * * Author: Verena Kaynig * * Kontakt: verena.kaynig@inf.ethz.ch * * * * SIFT implementation: Stephan Saalfeld * * **************************************************************************** */ package lenscorrection; import java.awt.Color; import java.awt.geom.AffineTransform; import java.io.BufferedWriter; import java.io.File; import java.io.FileWriter; import java.io.FilenameFilter; import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; import java.util.Collections; import java.util.List; import java.util.concurrent.atomic.AtomicInteger; import org.jfree.chart.ChartFactory; import org.jfree.chart.JFreeChart; import org.jfree.chart.plot.PlotOrientation; import org.jfree.data.category.DefaultCategoryDataset; import Jama.Matrix; import ij.IJ; import ij.ImagePlus; import ij.gui.GenericDialog; import ij.io.DirectoryChooser; import ij.io.FileSaver; import ij.io.Opener; import ij.plugin.PlugIn; import ij.process.ColorProcessor; import ij.process.ImageProcessor; import mpi.fruitfly.general.MultiThreading; import mpicbg.ij.SIFT; import mpicbg.imagefeatures.Feature; import mpicbg.imagefeatures.FloatArray2DSIFT; import mpicbg.models.AbstractAffineModel2D; import mpicbg.models.AffineModel2D; import mpicbg.models.NoninvertibleModelException; import mpicbg.models.NotEnoughDataPointsException; import mpicbg.models.PointMatch; import mpicbg.models.RigidModel2D; import mpicbg.models.SimilarityModel2D; import mpicbg.models.TranslationModel2D; public class Distortion_Correction implements PlugIn { static public class PointMatchCollectionAndAffine { final public AffineTransform affine; final public Collection<PointMatch> pointMatches; public PointMatchCollectionAndAffine(final AffineTransform affine, final Collection<PointMatch> pointMatches) { this.affine = affine; this.pointMatches = pointMatches; } } static public class BasicParam { final public FloatArray2DSIFT.Param sift = new FloatArray2DSIFT.Param(); /** * Closest/next closest neighbour distance ratio */ public float rod = 0.92f; /** * Maximal allowed alignment error in px */ public float maxEpsilon = 32.0f; /** * Inlier/candidates ratio */ public float minInlierRatio = 0.2f; /** * Implemeted transformation models for choice */ final static public String[] modelStrings = new String[] { "Translation", "Rigid", "Similarity", "Affine" }; public int expectedModelIndex = 1; /** * Order of the polynomial kernel */ public int dimension = 5; /** * Regularization factor */ public double lambda = 0.01; public BasicParam() { sift.fdSize = 8; sift.maxOctaveSize = 1600; sift.minOctaveSize = 400; } public void addFields(final GenericDialog gd) { SIFT.addFields(gd, sift); gd.addNumericField("closest/next_closest_ratio :", rod, 2); gd.addMessage("Geometric Consensus Filter :"); gd.addNumericField("maximal_alignment_error :", maxEpsilon, 2, 6, "px"); gd.addNumericField("inlier_ratio :", minInlierRatio, 2); gd.addChoice("expected_transformation :", modelStrings, modelStrings[expectedModelIndex]); gd.addMessage("Lens Model :"); gd.addNumericField("power_of_polynomial_kernel :", dimension, 0); gd.addNumericField("lambda :", lambda, 6); } public boolean readFields(final GenericDialog gd) { SIFT.readFields(gd, sift); rod = (float) gd.getNextNumber(); maxEpsilon = (float) gd.getNextNumber(); minInlierRatio = (float) gd.getNextNumber(); expectedModelIndex = gd.getNextChoiceIndex(); dimension = (int) gd.getNextNumber(); lambda = (double) gd.getNextNumber(); return !gd.invalidNumber(); } public boolean setup(final String title) { final GenericDialog gd = new GenericDialog(title); addFields(gd); do { gd.showDialog(); if (gd.wasCanceled()) return false; } while (!readFields(gd)); return true; } @Override public BasicParam clone() { final BasicParam p = new BasicParam(); p.sift.set(sift); p.dimension = dimension; p.expectedModelIndex = expectedModelIndex; p.lambda = lambda; p.maxEpsilon = maxEpsilon; p.minInlierRatio = minInlierRatio; p.rod = rod; return p; } } static public class PluginParam extends BasicParam { public String source_dir = ""; public String target_dir = ""; public String saveFileName = "distCorr.txt"; public boolean applyCorrection = true; public boolean visualizeResults = true; public int numberOfImages = 9; public int firstImageIndex = 0; /** * Original and calibrated images are supposed to have the same names, * but are in different directories */ public String[] names; public int saveOrLoad = 0; @Override public void addFields(final GenericDialog gd) { super.addFields(gd); } @Override public boolean readFields(final GenericDialog gd) { super.readFields(gd); return !gd.invalidNumber(); } /** * Setup as a three step dialog. */ @Override public boolean setup(final String title) { source_dir = ""; while (source_dir == "") { final DirectoryChooser dc = new DirectoryChooser("Calibration Images"); source_dir = dc.getDirectory(); if (null == source_dir) return false; source_dir = source_dir.replace('\\', '/'); if (!source_dir.endsWith("/")) source_dir += "/"; } final String exts = ".tif.jpg.png.gif.tiff.jpeg.bmp.pgm"; names = new File(source_dir).list(new FilenameFilter() { @Override public boolean accept(final File dir, final String name) { final int idot = name.lastIndexOf('.'); if (-1 == idot) return false; return exts.contains(name.substring(idot).toLowerCase()); } }); Arrays.sort(names); final GenericDialog gd = new GenericDialog(title); gd.addNumericField("number_of_images :", 9, 0); gd.addChoice("first_image :", names, names[0]); gd.addNumericField("power_of_polynomial_kernel :", dimension, 0); gd.addNumericField("lambda :", lambda, 6); gd.addCheckbox("apply_correction_to_images", applyCorrection); gd.addCheckbox("visualize results", visualizeResults); final String[] options = new String[] { "save", "load" }; gd.addChoice("What to do? ", options, options[saveOrLoad]); gd.addStringField("file_name: ", saveFileName); gd.showDialog(); if (gd.wasCanceled()) return false; numberOfImages = (int) gd.getNextNumber(); firstImageIndex = gd.getNextChoiceIndex(); dimension = (int) gd.getNextNumber(); lambda = gd.getNextNumber(); applyCorrection = gd.getNextBoolean(); visualizeResults = gd.getNextBoolean(); saveOrLoad = gd.getNextChoiceIndex(); saveFileName = gd.getNextString(); if (saveOrLoad == 0 || visualizeResults) { final GenericDialog gds = new GenericDialog(title); SIFT.addFields(gds, sift); gds.addNumericField("closest/next_closest_ratio :", rod, 2); gds.addMessage("Geometric Consensus Filter:"); gds.addNumericField("maximal_alignment_error :", maxEpsilon, 2, 6, "px"); gds.addNumericField("inlier_ratio :", minInlierRatio, 2); gds.addChoice("expected_transformation :", modelStrings, modelStrings[expectedModelIndex]); gds.showDialog(); if (gds.wasCanceled()) return false; SIFT.readFields(gds, sift); rod = (float) gds.getNextNumber(); maxEpsilon = (float) gds.getNextNumber(); minInlierRatio = (float) gds.getNextNumber(); expectedModelIndex = gds.getNextChoiceIndex(); return !(gd.invalidNumber() || gds.invalidNumber()); } return !gd.invalidNumber(); } } static final public BasicParam p = new BasicParam(); static final public PluginParam sp = new PluginParam(); NonLinearTransform nlt = new NonLinearTransform(); AbstractAffineModel2D<?>[] models; @Override public void run(final String arg) { if (!sp.setup("Lens Correction")) return; List<List<PointMatch>> inliers = null; if (sp.saveOrLoad == 0) { //IJ.log( sp.source_dir + sp.names[ 0 ] ); final ImagePlus imgTmp = new Opener().openImage(sp.source_dir + sp.names[0]); final int imageWidth = imgTmp.getWidth(), imageHeight = imgTmp.getHeight(); /** imgTmp was just needed to get width and height of the images */ imgTmp.flush(); final List<Feature>[] siftFeatures = extractSIFTFeaturesThreaded(sp.numberOfImages, sp.source_dir, sp.names); final List<PointMatch>[] inliersTmp = new ArrayList[sp.numberOfImages * (sp.numberOfImages - 1)]; models = new AbstractAffineModel2D[sp.numberOfImages * (sp.numberOfImages - 1)]; IJ.showStatus("Estimating Correspondences"); for (int i = 0; i < sp.numberOfImages; ++i) { IJ.log("Estimating correspondences of image " + i); IJ.showProgress((i + 1), sp.numberOfImages); extractSIFTPointsThreaded(i, siftFeatures, inliersTmp, models); } int wholeCount = 0; inliers = new ArrayList<List<PointMatch>>(); for (int i = 0; i < inliersTmp.length; i++) { // if vector at inliersTmp[i] contains only one null element, // its size is still 1 if (inliersTmp[i].size() > 10) { wholeCount += inliersTmp[i].size(); } //if I do not do this then the models have not the //right index corresponding to the inliers positions in the vector inliers.add(inliersTmp[i]); } //this is really really ugly final double[][] tp = new double[wholeCount][6]; final double h1[][] = new double[wholeCount][2]; final double h2[][] = new double[wholeCount][2]; int count = 0; for (int i = 0; i < inliers.size(); ++i) { if (inliers.get(i).size() > 10) { final double[][] points1 = new double[inliers.get(i).size()][2]; final double[][] points2 = new double[inliers.get(i).size()][2]; for (int j = 0; j < inliers.get(i).size(); j++) { final double[] tmp1 = ((PointMatch) inliers.get(i).get(j)).getP1().getL(); final double[] tmp2 = ((PointMatch) inliers.get(i).get(j)).getP2().getL(); points1[j][0] = (double) tmp1[0]; points1[j][1] = (double) tmp1[1]; points2[j][0] = (double) tmp2[0]; points2[j][1] = (double) tmp2[1]; h1[count] = new double[] { (double) tmp1[0], (double) tmp1[1] }; h2[count] = new double[] { (double) tmp2[0], (double) tmp2[1] }; models[i].createAffine().getMatrix(tp[count]); count++; } } } //if ( sp.saveOrLoad == 0 ) //{ nlt = distortionCorrection(h1, h2, tp, sp.dimension, sp.lambda, imageWidth, imageHeight); nlt.visualizeSmall(sp.lambda); while (true) { final GenericDialog gdl = new GenericDialog("New lambda?"); gdl.addMessage( "If the distortion field shows a clear translation, \n it is likely that you need to increase lambda."); gdl.addNumericField("lambda :", sp.lambda, 6); gdl.showDialog(); if (gdl.wasCanceled()) break; sp.lambda = gdl.getNextNumber(); nlt = distortionCorrection(h1, h2, tp, sp.dimension, sp.lambda, imageWidth, imageHeight); nlt.visualizeSmall(sp.lambda); } nlt.save(sp.source_dir + sp.saveFileName); } //after all preprocessing is done, estimate the distortion correction transform if (sp.saveOrLoad == 1) { nlt.load(sp.source_dir + sp.saveFileName); nlt.print(); } if (sp.applyCorrection || sp.visualizeResults) { sp.target_dir = correctImages(); } if (sp.visualizeResults && (sp.saveOrLoad == 0)) { IJ.log("call nlt.visualize()"); nlt.visualize(); if (null != sp.target_dir) { IJ.log("Evaluating distortion correction..."); evaluateCorrection(inliers); } } IJ.log(" Done "); } public void evaluateCorrection(final List<List<PointMatch>> inliers) { IJ.showStatus("Evaluating Distortion Correction"); final double[][] original = new double[sp.numberOfImages][2]; final double[][] corrected = new double[sp.numberOfImages][2]; for (int i = sp.firstImageIndex; i < sp.numberOfImages; i++) { original[i] = evaluateCorrectionXcorr(i, sp.source_dir); corrected[i] = evaluateCorrectionXcorr(i, sp.target_dir); } final DefaultCategoryDataset dataset = new DefaultCategoryDataset(); final DefaultCategoryDataset datasetGain = new DefaultCategoryDataset(); final DefaultCategoryDataset datasetGrad = new DefaultCategoryDataset(); for (int i = 0; i < (original.length); ++i) { dataset.setValue(Math.abs(original[i][0]), "before", "image" + i); dataset.setValue(Math.abs(corrected[i][0]), "after", "image" + i); datasetGrad.setValue(Math.abs(original[i][1]), "before", "image" + i); datasetGrad.setValue(Math.abs(corrected[i][1]), "after", "image" + i); datasetGain.setValue(Math.abs(corrected[i][0]) - Math.abs(original[i][0]), "gray", "image" + i); datasetGain.setValue(Math.abs(corrected[i][1]) - Math.abs(original[i][1]), "grad", "image" + i); } final JFreeChart chart = ChartFactory.createBarChart("Xcorr before and after correction", "ImageNumber", "Xcorr", dataset, PlotOrientation.VERTICAL, false, true, false); final ImagePlus imp = new ImagePlus("Xcorr before and after correction Plot", chart.createBufferedImage(500, 300)); imp.show(); final JFreeChart chartGrad = ChartFactory.createBarChart("XcorrGradient before and after correction", "ImageNumber", "Xcorr", datasetGrad, PlotOrientation.VERTICAL, false, true, false); final ImagePlus impGrad = new ImagePlus("XcorrGradient before and after correction Plot", chartGrad.createBufferedImage(500, 300)); impGrad.show(); final JFreeChart chartGain = ChartFactory.createBarChart("Gain in Xcorr", "ImageNumber", "Xcorr", datasetGain, PlotOrientation.VERTICAL, false, true, false); final ImagePlus impGain = new ImagePlus("Gain in Xcorr Plot", chartGain.createBufferedImage(500, 300)); impGain.show(); visualizePoints(inliers); //write xcorr data to file String original0 = "", original1 = "", corrected0 = "", corrected1 = "", gain0 = "", gain1 = ""; for (int i = 0; i < (original.length); ++i) { original0 = original0 + Double.toString(original[i][0]) + "; "; original1 = original1 + Double.toString(original[i][1]) + "; "; corrected0 = corrected0 + Double.toString(corrected[i][0]) + "; "; corrected1 = corrected1 + Double.toString(corrected[i][1]) + "; "; gain0 = gain0 + Double.toString(Math.abs(corrected[i][0]) - Math.abs(original[i][0])) + "; "; gain1 = gain1 + Double.toString(Math.abs(corrected[i][1]) - Math.abs(original[i][1])) + "; "; } try { final BufferedWriter out = new BufferedWriter(new FileWriter(sp.source_dir + "xcorrData.log")); out.write(original0); out.newLine(); out.newLine(); out.write(original1); out.newLine(); out.newLine(); out.write(corrected0); out.newLine(); out.newLine(); out.write(corrected1); out.newLine(); out.newLine(); out.write(gain0); out.newLine(); out.newLine(); out.write(gain1); out.newLine(); out.close(); } catch (final Exception e) { System.err.println("Error: " + e.getMessage()); } ; } protected void extractSIFTPoints(final int index, final List<Feature>[] siftFeatures, final List<List<PointMatch>> inliers, final List<AbstractAffineModel2D<?>> models) { //save all matching candidates final List<List<PointMatch>> candidates = new ArrayList<List<PointMatch>>(); for (int j = 0; j < siftFeatures.length; j++) { if (index == j) continue; candidates.add(FloatArray2DSIFT.createMatches(siftFeatures[index], siftFeatures[j], 1.5f, null, Float.MAX_VALUE, 0.5f)); } //get rid of the outliers and save the transformations to match the inliers for (int i = 0; i < candidates.size(); ++i) { final List<PointMatch> tmpInliers = new ArrayList<PointMatch>(); final AbstractAffineModel2D<?> m; switch (sp.expectedModelIndex) { case 0: m = new TranslationModel2D(); break; case 1: m = new RigidModel2D(); break; case 2: m = new SimilarityModel2D(); break; case 3: m = new AffineModel2D(); break; default: return; } try { m.filterRansac(candidates.get(i), tmpInliers, 1000, sp.maxEpsilon, sp.minInlierRatio, 10); } catch (final NotEnoughDataPointsException e) { e.printStackTrace(); } inliers.add(tmpInliers); models.add(m); } } /** * Estimates a polynomial distortion model for a set of * {@linkplain PointMatch corresponding points} and returns the inverse of * this model which then may be used to undo the distortion. * * @param pointMatches * a list of matches, plus affines that transfer each match * into a common image frame * @param dimension the order of the polynomial model * @param lambda regularization factor * @param imageWidth * @param imageHeight * * @return a polynomial distortion model to undo the distortion */ static public NonLinearTransform createInverseDistortionModel( final Collection<PointMatchCollectionAndAffine> pointMatches, final int dimension, final double lambda, final int imageWidth, final int imageHeight) { int wholeCount = 0; for (final PointMatchCollectionAndAffine pma : pointMatches) if (pma.pointMatches.size() > 10) wholeCount += pma.pointMatches.size(); final double[][] tp = new double[wholeCount][6]; final double h1[][] = new double[wholeCount][2]; final double h2[][] = new double[wholeCount][2]; int count = 0; for (final PointMatchCollectionAndAffine pma : pointMatches) { if (pma.pointMatches.size() > 10) { int i = 0; for (final PointMatch match : pma.pointMatches) { final double[] tmp1 = match.getP1().getL(); final double[] tmp2 = match.getP2().getL(); h1[count] = new double[] { (double) tmp1[0], (double) tmp1[1] }; h2[count] = new double[] { (double) tmp2[0], (double) tmp2[1] }; pma.affine.getMatrix(tp[count]); count++; ++i; } } } final NonLinearTransform nlt = distortionCorrection(h1, h2, tp, dimension, lambda, imageWidth, imageHeight); nlt.inverseTransform(h1); return nlt; } static protected NonLinearTransform distortionCorrection(final double hack1[][], final double hack2[][], final double transformParams[][], final int dimension, final double lambda, final int w, final int h) { IJ.showStatus("Getting the Distortion Field"); final NonLinearTransform nlt = new NonLinearTransform(dimension, w, h); nlt.estimateDistortion(hack1, hack2, transformParams, lambda, w, h); nlt.inverseTransform(hack1); return nlt; } protected String correctImages() { if (!sp.applyCorrection) { sp.target_dir = System.getProperty("user.dir").replace('\\', '/') + "/distCorr_tmp/"; System.out.println("Tmp target directory: " + sp.target_dir); if (new File(sp.target_dir).exists()) { System.out.println("removing old tmp directory!"); final String[] filesToDelete = new File(sp.target_dir).list(); for (int i = 0; i < filesToDelete.length; i++) { System.out.println(filesToDelete[i]); final boolean deleted = new File(sp.target_dir + filesToDelete[i]).delete(); if (!deleted) IJ.log("Error: Could not remove temporary directory!"); } new File(sp.target_dir).delete(); } try { // Create one directory final boolean success = (new File(sp.target_dir)).mkdir(); if (success) new File(sp.target_dir).deleteOnExit(); } catch (final Exception e) { IJ.showMessage("Error! Could not create temporary directory. " + e.getMessage()); } } if (sp.target_dir == "" || null == sp.target_dir) { final DirectoryChooser dc = new DirectoryChooser("Target Directory"); sp.target_dir = dc.getDirectory(); if (null == sp.target_dir) return null; sp.target_dir = sp.target_dir.replace('\\', '/'); if (!sp.target_dir.endsWith("/")) sp.target_dir += "/"; } final String[] namesTarget = new File(sp.target_dir).list(new FilenameFilter() { @Override public boolean accept(final File dir, final String namesTarget) { final int idot = namesTarget.lastIndexOf('.'); if (-1 == idot) return false; return namesTarget.contains(namesTarget.substring(idot).toLowerCase()); } }); if (namesTarget.length > 0) IJ.showMessage("Overwrite Message", "There are already images in that directory. These will be used for evaluation."); else { IJ.showStatus("Correcting Images"); final Thread[] threads = MultiThreading.newThreads(); final AtomicInteger ai = new AtomicInteger(sp.applyCorrection ? 0 : sp.firstImageIndex); for (int ithread = 0; ithread < threads.length; ++ithread) { threads[ithread] = new Thread() { @Override public void run() { setPriority(Thread.NORM_PRIORITY); for (int i = ai.getAndIncrement(); i < (sp.applyCorrection ? sp.names.length : (sp.firstImageIndex + sp.numberOfImages)); i = ai.getAndIncrement()) { IJ.log("Correcting image " + sp.names[i]); final ImagePlus imps = new Opener().openImage(sp.source_dir + sp.names[i]); imps.setProcessor(imps.getTitle(), imps.getProcessor().convertToShort(false)); final ImageProcessor[] transErg = nlt.transform(imps.getProcessor()); imps.setProcessor(imps.getTitle(), transErg[0]); if (!sp.applyCorrection) new File(sp.target_dir + sp.names[i]).deleteOnExit(); new FileSaver(imps).saveAsTiff(sp.target_dir + sp.names[i]); } } }; } MultiThreading.startAndJoin(threads); } return sp.target_dir; } protected double[] evaluateCorrectionXcorr(final int index, final String directory) { final ImagePlus im1 = new Opener().openImage(directory + sp.names[index]); im1.setProcessor(sp.names[index], im1.getProcessor().convertToShort(false)); int count = 0; final ArrayList<Double> xcorrVals = new ArrayList<Double>(); final ArrayList<Double> xcorrValsGrad = new ArrayList<Double>(); for (int i = 0; i < sp.numberOfImages; i++) { if (i == index) { continue; } if (models[index * (sp.numberOfImages - 1) + count] == null) { count++; continue; } final ImagePlus newImg = new Opener().openImage(directory + sp.names[i + sp.firstImageIndex]); newImg.setProcessor(newImg.getTitle(), newImg.getProcessor().convertToShort(false)); newImg.setProcessor(sp.names[i + sp.firstImageIndex], applyTransformToImageInverse( models[index * (sp.numberOfImages - 1) + count], newImg.getProcessor())); // If you want to see the stitching improvement run this // ImageProcessor testIp = im1.getProcessor().duplicate(); // for ( int x=0; x < testIp.getWidth(); x++){ // for (int y=0; y < testIp.getHeight(); y++){ // testIp.set(x, y, Math.abs(im1.getProcessor().get(x,y) - // newImg.getProcessor().get(x,y))); // } // } // ImagePlus testImg = new ImagePlus(sp.names[index] + " minus " + // sp.names[i], testIp); // testImg.show(); // im1.show(); // newImg.show(); xcorrVals.add(getXcorrBlackOut(im1.getProcessor(), newImg.getProcessor())); xcorrValsGrad.add(getXcorrBlackOutGradient(im1.getProcessor(), newImg.getProcessor())); count++; } Collections.sort(xcorrVals); Collections.sort(xcorrValsGrad); // double[] medians = { xcorrVals.get( xcorrVals.size() / 2 ), xcorrValsGrad.get( xcorrValsGrad.size() / 2 ) }; double m1 = 0.0, m2 = 0.0; for (int i = 0; i < xcorrVals.size(); i++) { m1 += xcorrVals.get(i); m2 += xcorrValsGrad.get(i); } m1 /= xcorrVals.size(); m2 /= xcorrVals.size(); final double[] means = { m1, m2 }; return means; //return medians; } ImageProcessor applyTransformToImageInverse(final AbstractAffineModel2D<?> a, final ImageProcessor ip) { final ImageProcessor newIp = ip.duplicate(); newIp.max(0.0); for (int x = 0; x < ip.getWidth(); x++) { for (int y = 0; y < ip.getHeight(); y++) { final double[] position = new double[] { x, y }; // float[] newPosition = a.apply(position); double[] newPosition = new double[] { 0, 0 }; try { newPosition = a.applyInverse(position); } catch (final NoninvertibleModelException e) { } final int xn = (int) newPosition[0]; final int yn = (int) newPosition[1]; if ((xn >= 0) && (yn >= 0) && (xn < ip.getWidth()) && (yn < ip.getHeight())) newIp.set(xn, yn, ip.get(x, y)); } } return newIp; } double getXcorrBlackOutGradient(final ImageProcessor ip1, final ImageProcessor ip2) { final ImageProcessor ip1g = getGradientSobel(ip1); final ImageProcessor ip2g = getGradientSobel(ip2); return getXcorrBlackOut(ip1g, ip2g); } //this blends out gradients that include black pixels to make the sharp border caused //by the nonlinear transformation not disturb the gradient comparison //FIXME: this should be handled by a mask image! ImageProcessor getGradientSobel(final ImageProcessor ip) { final ImageProcessor ipGrad = ip.duplicate(); ipGrad.max(0.0); for (int i = 1; i < ipGrad.getWidth() - 1; i++) { for (int j = 1; j < ipGrad.getHeight() - 1; j++) { if (ip.get(i - 1, j - 1) == 0 || ip.get(i - 1, j) == 0 || ip.get(i - 1, j + 1) == 0 || ip.get(i, j - 1) == 0 || ip.get(i, j) == 0 || ip.get(i, j + 1) == 0 || ip.get(i + 1, j - 1) == 0 || ip.get(i + 1, j) == 0 || ip.get(i + 1, j + 1) == 0) continue; final double gradX = (double) -ip.get(i - 1, j - 1) - 2 * ip.get(i - 1, j) - ip.get(i - 1, j + 1) + ip.get(i + 1, j - 1) + 2 * ip.get(i + 1, j) + ip.get(i + 1, j + 1); final double gradY = (double) -ip.get(i - 1, j - 1) - 2 * ip.get(i, j - 1) - ip.get(i + 1, j - 1) + ip.get(i - 1, j + 1) + 2 * ip.get(i, j + 1) + ip.get(i + 1, j + 1); final double mag = Math.sqrt(gradX * gradX + gradY * gradY); ipGrad.setf(i, j, (float) mag); } } return ipGrad; } //tested the result against matlab routine, this worked fine static double getXcorrBlackOut(ImageProcessor ip1, ImageProcessor ip2) { ip1 = ip1.convertToFloat(); ip2 = ip2.convertToFloat(); //If this is not done, the black area from the transformed image influences xcorr result //better alternative would be to use mask images and only calculate xcorr of //the region present in both images. for (int i = 0; i < ip1.getWidth(); i++) { for (int j = 0; j < ip1.getHeight(); j++) { if (ip1.get(i, j) == 0) ip2.set(i, j, 0); if (ip2.get(i, j) == 0) ip1.set(i, j, 0); } } // FloatProcessor ip1f = (FloatProcessor)ip1.convertToFloat(); // FloatProcessor ip2f = (FloatProcessor)ip2.convertToFloat(); final float[] data1 = (float[]) ip1.getPixels(); final float[] data2 = (float[]) ip2.getPixels(); final double[] data1b = new double[data1.length]; final double[] data2b = new double[data2.length]; int count = 0; double mean1 = 0.0, mean2 = 0.0; for (int i = 0; i < data1.length; i++) { //if ((data1[i] == 0) || (data2[i] == 0)) //continue; data1b[i] = data1[i]; data2b[i] = data2[i]; mean1 += data1b[i]; mean2 += data2b[i]; count++; } mean1 /= (double) count; mean2 /= (double) count; double L2_1 = 0.0, L2_2 = 0.0; for (int i = 0; i < count; i++) { L2_1 += (data1b[i] - mean1) * (data1b[i] - mean1); L2_2 += (data2b[i] - mean2) * (data2b[i] - mean2); } L2_1 = Math.sqrt(L2_1); L2_2 = Math.sqrt(L2_2); double xcorr = 0.0; for (int i = 0; i < count; i++) { xcorr += ((data1b[i] - mean1) / L2_1) * ((data2b[i] - mean2) / L2_2); } //System.out.println("XcorrVal: " + xcorr); return xcorr; } void visualizePoints(final List<List<PointMatch>> inliers) { final ColorProcessor ip = new ColorProcessor(nlt.getWidth(), nlt.getHeight()); ip.setColor(Color.red); ip.setLineWidth(5); for (int i = 0; i < inliers.size(); i++) { for (int j = 0; j < inliers.get(i).size(); j++) { final double[] tmp1 = inliers.get(i).get(j).getP1().getW(); final double[] tmp2 = inliers.get(i).get(j).getP2().getL(); ip.setColor(Color.red); ip.drawDot((int) tmp2[0], (int) tmp2[1]); ip.setColor(Color.blue); ip.drawDot((int) tmp1[0], (int) tmp1[1]); } } final ImagePlus points = new ImagePlus("Corresponding Points after correction", ip); points.show(); } public void getTransform(final double[][] points1, final double[][] points2, final double[][] transformParams) { final double[][] p1 = new double[points1.length][3]; final double[][] p2 = new double[points2.length][3]; for (int i = 0; i < points1.length; i++) { p1[i][0] = points1[i][0]; p1[i][1] = points1[i][1]; p1[i][2] = 100.0; p2[i][0] = points2[i][0]; p2[i][1] = points2[i][1]; p2[i][2] = 100.0; } final Matrix s1 = new Matrix(p1); final Matrix s2 = new Matrix(p2); Matrix t = (s1.transpose().times(s1)).inverse().times(s1.transpose()).times(s2); t = t.inverse(); for (int i = 0; i < transformParams.length; i++) { if (transformParams[i][0] == -10) { transformParams[i][0] = t.get(0, 0); transformParams[i][1] = t.get(0, 1); transformParams[i][2] = t.get(1, 0); transformParams[i][3] = t.get(1, 1); transformParams[i][4] = t.get(2, 0); transformParams[i][5] = t.get(2, 1); } } t.print(1, 1); } static List<Feature>[] extractSIFTFeaturesThreaded(final int numberOfImages, final String directory, final String[] names) { //extract all SIFT Features final List<Feature>[] siftFeatures = new ArrayList[numberOfImages]; final Thread[] threads = MultiThreading.newThreads(); final AtomicInteger ai = new AtomicInteger(0); // start at second slice IJ.showStatus("Extracting SIFT Features"); for (int ithread = 0; ithread < threads.length; ++ithread) { threads[ithread] = new Thread() { @Override public void run() { for (int i = ai.getAndIncrement(); i < numberOfImages; i = ai.getAndIncrement()) { final ArrayList<Feature> fs = new ArrayList<Feature>(); final ImagePlus imps = new Opener().openImage(directory + names[i + sp.firstImageIndex]); imps.setProcessor(imps.getTitle(), imps.getProcessor().convertToFloat()); final FloatArray2DSIFT sift = new FloatArray2DSIFT(sp.sift.clone()); final SIFT ijSIFT = new SIFT(sift); ijSIFT.extractFeatures(imps.getProcessor(), fs); Collections.sort(fs); IJ.log("Extracting SIFT of image: " + i); siftFeatures[i] = fs; } } }; } MultiThreading.startAndJoin(threads); return siftFeatures; } static protected void extractSIFTPointsThreaded(final int index, final List<Feature>[] siftFeatures, final List<PointMatch>[] inliers, final AbstractAffineModel2D<?>[] models) { // save all matching candidates final List<PointMatch>[] candidates = new List[siftFeatures.length - 1]; final Thread[] threads = MultiThreading.newThreads(); final AtomicInteger ai = new AtomicInteger(0); // start at second // slice for (int ithread = 0; ithread < threads.length; ++ithread) { threads[ithread] = new Thread() { @Override public void run() { setPriority(Thread.NORM_PRIORITY); for (int j = ai.getAndIncrement(); j < candidates.length; j = ai.getAndIncrement()) { final int i = (j < index ? j : j + 1); candidates[j] = FloatArray2DSIFT.createMatches(siftFeatures[index], siftFeatures[i], 1.5f, null, Float.MAX_VALUE, 0.5f); } } }; } MultiThreading.startAndJoin(threads); // get rid of the outliers and save the rigid transformations to match // the inliers final AtomicInteger ai2 = new AtomicInteger(0); for (int ithread = 0; ithread < threads.length; ++ithread) { threads[ithread] = new Thread() { @Override public void run() { setPriority(Thread.NORM_PRIORITY); for (int i = ai2.getAndIncrement(); i < candidates.length; i = ai2.getAndIncrement()) { final List<PointMatch> tmpInliers = new ArrayList<PointMatch>(); // RigidModel2D m = // RigidModel2D.estimateBestModel(candidates.get(i), // tmpInliers, sp.min_epsilon, sp.max_epsilon, // sp.min_inlier_ratio); final AbstractAffineModel2D<?> m; switch (sp.expectedModelIndex) { case 0: m = new TranslationModel2D(); break; case 1: m = new RigidModel2D(); break; case 2: m = new SimilarityModel2D(); break; case 3: m = new AffineModel2D(); break; default: return; } boolean modelFound = false; try { modelFound = m.filterRansac(candidates[i], tmpInliers, 1000, sp.maxEpsilon, sp.minInlierRatio, 10); } catch (final NotEnoughDataPointsException e) { modelFound = false; } if (modelFound) IJ.log("Model found:\n " + candidates[i].size() + " candidates\n " + tmpInliers.size() + " inliers\n " + String.format("%.2f", m.getCost()) + "px average displacement"); else IJ.log("No Model found."); inliers[index * (sp.numberOfImages - 1) + i] = tmpInliers; models[index * (sp.numberOfImages - 1) + i] = m; // System.out.println("**** MODEL ADDED: " + // (index*(sp.numberOfImages-1)+i)); } } }; } MultiThreading.startAndJoin(threads); } }