Search in sources :

Example 1 with MovingLeastSquaresTransform2

use of mpicbg.trakem2.transform.MovingLeastSquaresTransform2 in project TrakEM2 by trakem2.

the class ElasticLayerAlignment method exec.

/**
 * @param param
 * @param layerRange
 * @param fixedLayers
 * @param emptyLayers
 * @param box
 * @param filter
 * @param useTps true if using TPS transforms, otherwise MLS
 * @throws Exception
 */
@SuppressWarnings("deprecation")
public final void exec(final Param param, final Project project, final List<Layer> layerRange, final Set<Layer> fixedLayers, final Set<Layer> emptyLayers, final Rectangle box, final boolean propagateTransformBefore, final boolean propagateTransformAfter, final Filter<Patch> filter) throws Exception {
    final ExecutorService service = ExecutorProvider.getExecutorService(1.0f);
    /* create tiles and models for all layers */
    final ArrayList<Tile<?>> tiles = new ArrayList<Tile<?>>();
    for (int i = 0; i < layerRange.size(); ++i) {
        switch(param.desiredModelIndex) {
            case 0:
                tiles.add(new Tile<TranslationModel2D>(new TranslationModel2D()));
                break;
            case 1:
                tiles.add(new Tile<RigidModel2D>(new RigidModel2D()));
                break;
            case 2:
                tiles.add(new Tile<SimilarityModel2D>(new SimilarityModel2D()));
                break;
            case 3:
                tiles.add(new Tile<AffineModel2D>(new AffineModel2D()));
                break;
            case 4:
                tiles.add(new Tile<HomographyModel2D>(new HomographyModel2D()));
                break;
            default:
                return;
        }
    }
    /* collect all pairs of slices for which a model could be found */
    final ArrayList<Triple<Integer, Integer, AbstractModel<?>>> pairs = new ArrayList<Triple<Integer, Integer, AbstractModel<?>>>();
    if (!param.isAligned) {
        preAlignStack(param, project, layerRange, box, filter, pairs);
    } else {
        for (int i = 0; i < layerRange.size(); ++i) {
            final int range = Math.min(layerRange.size(), i + param.maxNumNeighbors + 1);
            for (int j = i + 1; j < range; ++j) {
                pairs.add(new Triple<Integer, Integer, AbstractModel<?>>(i, j, new TranslationModel2D()));
            }
        }
    }
    /* Elastic alignment */
    /* Initialization */
    final TileConfiguration initMeshes = new TileConfiguration();
    final int meshWidth = (int) Math.ceil(box.width * param.layerScale);
    final int meshHeight = (int) Math.ceil(box.height * param.layerScale);
    final ArrayList<SpringMesh> meshes = new ArrayList<SpringMesh>(layerRange.size());
    for (int i = 0; i < layerRange.size(); ++i) {
        meshes.add(new SpringMesh(param.resolutionSpringMesh, meshWidth, meshHeight, param.stiffnessSpringMesh, param.maxStretchSpringMesh * param.layerScale, param.dampSpringMesh));
    }
    // final int blockRadius = Math.max( 32, meshWidth / p.resolutionSpringMesh / 2 );
    final int blockRadius = Math.max(16, mpicbg.util.Util.roundPos(param.layerScale * param.blockRadius));
    Utils.log("effective block radius = " + blockRadius);
    final ArrayList<Future<BlockMatchPairCallable.BlockMatchResults>> futures = new ArrayList<Future<BlockMatchPairCallable.BlockMatchResults>>(pairs.size());
    for (final Triple<Integer, Integer, AbstractModel<?>> pair : pairs) {
        /* free memory */
        project.getLoader().releaseAll();
        final SpringMesh m1 = meshes.get(pair.a);
        final SpringMesh m2 = meshes.get(pair.b);
        final ArrayList<Vertex> v1 = m1.getVertices();
        final ArrayList<Vertex> v2 = m2.getVertices();
        final Layer layer1 = layerRange.get(pair.a);
        final Layer layer2 = layerRange.get(pair.b);
        final boolean layer1Fixed = fixedLayers.contains(layer1);
        final boolean layer2Fixed = fixedLayers.contains(layer2);
        if (!(layer1Fixed && layer2Fixed)) {
            final BlockMatchPairCallable bmpc = new BlockMatchPairCallable(pair, layerRange, layer1Fixed, layer2Fixed, filter, param, v1, v2, box);
            futures.add(service.submit(bmpc));
        }
    }
    for (final Future<BlockMatchPairCallable.BlockMatchResults> future : futures) {
        final BlockMatchPairCallable.BlockMatchResults results = future.get();
        final Collection<PointMatch> pm12 = results.pm12, pm21 = results.pm21;
        final Triple<Integer, Integer, AbstractModel<?>> pair = results.pair;
        final Tile<?> t1 = tiles.get(pair.a);
        final Tile<?> t2 = tiles.get(pair.b);
        final SpringMesh m1 = meshes.get(pair.a);
        final SpringMesh m2 = meshes.get(pair.b);
        final double springConstant = 1.0 / (pair.b - pair.a);
        final boolean layer1Fixed = results.layer1Fixed;
        final boolean layer2Fixed = results.layer2Fixed;
        if (layer1Fixed) {
            initMeshes.fixTile(t1);
        } else {
            if (param.useLocalSmoothnessFilter) {
                Utils.log(pair.a + " > " + pair.b + ": " + pm12.size() + " candidates passed local smoothness filter.");
            } else {
                Utils.log(pair.a + " > " + pair.b + ": found " + pm12.size() + " correspondences.");
            }
            for (final PointMatch pm : pm12) {
                final Vertex p1 = (Vertex) pm.getP1();
                final Vertex p2 = new Vertex(pm.getP2());
                p1.addSpring(p2, new Spring(0, springConstant));
                m2.addPassiveVertex(p2);
            }
            /*
                * adding Tiles to the initialing TileConfiguration, adding a Tile
                * multiple times does not harm because the TileConfiguration is
                * backed by a Set.
                */
            if (pm12.size() > pair.c.getMinNumMatches()) {
                initMeshes.addTile(t1);
                initMeshes.addTile(t2);
                t1.connect(t2, pm12);
            }
        }
        if (layer2Fixed)
            initMeshes.fixTile(t2);
        else {
            if (param.useLocalSmoothnessFilter) {
                Utils.log(pair.a + " < " + pair.b + ": " + pm21.size() + " candidates passed local smoothness filter.");
            } else {
                Utils.log(pair.a + " < " + pair.b + ": found " + pm21.size() + " correspondences.");
            }
            for (final PointMatch pm : pm21) {
                final Vertex p1 = (Vertex) pm.getP1();
                final Vertex p2 = new Vertex(pm.getP2());
                p1.addSpring(p2, new Spring(0, springConstant));
                m1.addPassiveVertex(p2);
            }
            /*
                * adding Tiles to the initialing TileConfiguration, adding a Tile
                * multiple times does not harm because the TileConfiguration is
                * backed by a Set.
                */
            if (pm21.size() > pair.c.getMinNumMatches()) {
                initMeshes.addTile(t1);
                initMeshes.addTile(t2);
                t2.connect(t1, pm21);
            }
        }
        Utils.log(pair.a + " <> " + pair.b + " spring constant = " + springConstant);
    }
    /* pre-align by optimizing a piecewise linear model */
    initMeshes.optimize(param.maxEpsilon * param.layerScale, param.maxIterationsSpringMesh, param.maxPlateauwidthSpringMesh);
    for (int i = 0; i < layerRange.size(); ++i) meshes.get(i).init(tiles.get(i).getModel());
    /* optimize the meshes */
    try {
        final long t0 = System.currentTimeMillis();
        Utils.log("Optimizing spring meshes...");
        if (param.useLegacyOptimizer) {
            Utils.log("  ...using legacy optimizer...");
            SpringMesh.optimizeMeshes2(meshes, param.maxEpsilon * param.layerScale, param.maxIterationsSpringMesh, param.maxPlateauwidthSpringMesh, param.visualize);
        } else {
            SpringMesh.optimizeMeshes(meshes, param.maxEpsilon * param.layerScale, param.maxIterationsSpringMesh, param.maxPlateauwidthSpringMesh, param.visualize);
        }
        Utils.log("Done optimizing spring meshes. Took " + (System.currentTimeMillis() - t0) + " ms");
    } catch (final NotEnoughDataPointsException e) {
        Utils.log("There were not enough data points to get the spring mesh optimizing.");
        e.printStackTrace();
        return;
    }
    /* translate relative to bounding box */
    for (final SpringMesh mesh : meshes) {
        for (final PointMatch pm : mesh.getVA().keySet()) {
            final Point p1 = pm.getP1();
            final Point p2 = pm.getP2();
            final double[] l = p1.getL();
            final double[] w = p2.getW();
            l[0] = l[0] / param.layerScale + box.x;
            l[1] = l[1] / param.layerScale + box.y;
            w[0] = w[0] / param.layerScale + box.x;
            w[1] = w[1] / param.layerScale + box.y;
        }
    }
    /* free memory */
    project.getLoader().releaseAll();
    final Layer first = layerRange.get(0);
    final List<Layer> layers = first.getParent().getLayers();
    final LayerSet ls = first.getParent();
    final Area infArea = AreaUtils.infiniteArea();
    final List<VectorData> vectorData = new ArrayList<VectorData>();
    for (final Layer layer : ls.getLayers()) {
        vectorData.addAll(Utils.castCollection(layer.getDisplayables(VectorData.class, false, true), VectorData.class, true));
    }
    vectorData.addAll(Utils.castCollection(ls.getZDisplayables(VectorData.class, true), VectorData.class, true));
    /* transfer layer transform into patch transforms and append to patches */
    if (propagateTransformBefore || propagateTransformAfter) {
        if (propagateTransformBefore) {
            final ThinPlateSplineTransform tps = makeTPS(meshes.get(0).getVA().keySet());
            final int firstLayerIndex = first.getParent().getLayerIndex(first.getId());
            for (int i = 0; i < firstLayerIndex; ++i) {
                applyTransformToLayer(layers.get(i), tps, filter);
                for (final VectorData vd : vectorData) {
                    vd.apply(layers.get(i), infArea, tps);
                }
            }
        }
        if (propagateTransformAfter) {
            final Layer last = layerRange.get(layerRange.size() - 1);
            final CoordinateTransform ct;
            if (param.useTps)
                ct = makeTPS(meshes.get(meshes.size() - 1).getVA().keySet());
            else {
                final MovingLeastSquaresTransform2 mls = new MovingLeastSquaresTransform2();
                mls.setMatches(meshes.get(meshes.size() - 1).getVA().keySet());
                ct = mls;
            }
            final int lastLayerIndex = last.getParent().getLayerIndex(last.getId());
            for (int i = lastLayerIndex + 1; i < layers.size(); ++i) {
                applyTransformToLayer(layers.get(i), ct, filter);
                for (final VectorData vd : vectorData) {
                    vd.apply(layers.get(i), infArea, ct);
                }
            }
        }
    }
    for (int l = 0; l < layerRange.size(); ++l) {
        IJ.showStatus("Applying transformation to patches ...");
        IJ.showProgress(0, layerRange.size());
        final Layer layer = layerRange.get(l);
        final ThinPlateSplineTransform tps = makeTPS(meshes.get(l).getVA().keySet());
        applyTransformToLayer(layer, tps, filter);
        for (final VectorData vd : vectorData) {
            vd.apply(layer, infArea, tps);
        }
        if (Thread.interrupted()) {
            Utils.log("Interrupted during applying transformations to patches.  No all patches have been updated.  Re-generate mipmaps manually.");
        }
        IJ.showProgress(l + 1, layerRange.size());
    }
    /* update patch mipmaps */
    final int firstLayerIndex;
    final int lastLayerIndex;
    if (propagateTransformBefore)
        firstLayerIndex = 0;
    else {
        firstLayerIndex = first.getParent().getLayerIndex(first.getId());
    }
    if (propagateTransformAfter)
        lastLayerIndex = layers.size() - 1;
    else {
        final Layer last = layerRange.get(layerRange.size() - 1);
        lastLayerIndex = last.getParent().getLayerIndex(last.getId());
    }
    for (int i = firstLayerIndex; i <= lastLayerIndex; ++i) {
        final Layer layer = layers.get(i);
        if (!(emptyLayers.contains(layer) || fixedLayers.contains(layer))) {
            for (final Patch patch : AlignmentUtils.filterPatches(layer, filter)) patch.updateMipMaps();
        }
    }
    Utils.log("Done.");
}
Also used : AbstractModel(mpicbg.models.AbstractModel) BlockMatchPairCallable(mpicbg.trakem2.align.concurrent.BlockMatchPairCallable) ArrayList(java.util.ArrayList) RigidModel2D(mpicbg.models.RigidModel2D) SimilarityModel2D(mpicbg.models.SimilarityModel2D) LayerSet(ini.trakem2.display.LayerSet) HomographyModel2D(mpicbg.models.HomographyModel2D) Spring(mpicbg.models.Spring) Triple(mpicbg.trakem2.util.Triple) Area(java.awt.geom.Area) TranslationModel2D(mpicbg.models.TranslationModel2D) TileConfiguration(mpicbg.models.TileConfiguration) CoordinateTransform(mpicbg.trakem2.transform.CoordinateTransform) Patch(ini.trakem2.display.Patch) NotEnoughDataPointsException(mpicbg.models.NotEnoughDataPointsException) Vertex(mpicbg.models.Vertex) SpringMesh(mpicbg.models.SpringMesh) ThinPlateSplineTransform(mpicbg.trakem2.transform.ThinPlateSplineTransform) VectorData(ini.trakem2.display.VectorData) AffineModel2D(mpicbg.models.AffineModel2D) Tile(mpicbg.models.Tile) Point(mpicbg.models.Point) Layer(ini.trakem2.display.Layer) Point(mpicbg.models.Point) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) PointMatch(mpicbg.models.PointMatch) MovingLeastSquaresTransform2(mpicbg.trakem2.transform.MovingLeastSquaresTransform2) ExecutorService(java.util.concurrent.ExecutorService) Future(java.util.concurrent.Future)

Example 2 with MovingLeastSquaresTransform2

use of mpicbg.trakem2.transform.MovingLeastSquaresTransform2 in project TrakEM2 by trakem2.

the class Align method createMLST.

/**
 * Temporary helper method that creates
 * @param matches
 * @param alpha
 * @return
 * @throws Exception
 */
public static final MovingLeastSquaresTransform2 createMLST(final Collection<PointMatch> matches, final double alpha) throws Exception {
    final MovingLeastSquaresTransform2 mlst = new MovingLeastSquaresTransform2();
    mlst.setAlpha(alpha);
    Class<? extends AbstractAffineModel2D<?>> c = AffineModel2D.class;
    switch(matches.size()) {
        case 1:
            c = TranslationModel2D.class;
            break;
        case 2:
            c = SimilarityModel2D.class;
            break;
        default:
            break;
    }
    mlst.setModel(c);
    mlst.setMatches(matches);
    return mlst;
}
Also used : MovingLeastSquaresTransform2(mpicbg.trakem2.transform.MovingLeastSquaresTransform2) AbstractAffineModel2D(mpicbg.models.AbstractAffineModel2D) AffineModel2D(mpicbg.models.AffineModel2D) InterpolatedAffineModel2D(mpicbg.models.InterpolatedAffineModel2D)

Example 3 with MovingLeastSquaresTransform2

use of mpicbg.trakem2.transform.MovingLeastSquaresTransform2 in project TrakEM2 by trakem2.

the class NonLinearTransformMode method createCT.

private CoordinateTransform createCT() throws Exception {
    final Collection<PointMatch> pm = new ArrayList<PointMatch>();
    for (final Point p : points) {
        pm.add(new PointMatch(new Point(p.getL()), new Point(p.getW())));
    }
    /*
		 * TODO replace this with the desired parameters of the transformation
		 */
    final MovingLeastSquaresTransform2 mlst = new MovingLeastSquaresTransform2();
    mlst.setAlpha(1.0f);
    Class<? extends AbstractAffineModel2D<?>> c = AffineModel2D.class;
    switch(points.size()) {
        case 1:
            c = TranslationModel2D.class;
            break;
        case 2:
            c = SimilarityModel2D.class;
            break;
        default:
            break;
    }
    mlst.setModel(c);
    mlst.setMatches(pm);
    return mlst;
}
Also used : PointMatch(mpicbg.models.PointMatch) MovingLeastSquaresTransform2(mpicbg.trakem2.transform.MovingLeastSquaresTransform2) ArrayList(java.util.ArrayList) AbstractAffineModel2D(mpicbg.models.AbstractAffineModel2D) AffineModel2D(mpicbg.trakem2.transform.AffineModel2D) Point(mpicbg.models.Point)

Example 4 with MovingLeastSquaresTransform2

use of mpicbg.trakem2.transform.MovingLeastSquaresTransform2 in project TrakEM2 by trakem2.

the class TestMovingLeastSquaresTransform2 method main.

public static final void main(String[] args) {
    try {
        // Preliminary test:
        System.out.println("Float.MIN_VALUE = " + Float.toString(Float.MIN_VALUE));
        System.out.println("Float.MAX_VALUE = " + Float.toString(Float.MAX_VALUE));
        System.out.println((float) Math.pow(10, -45));
        System.out.println(1.4f * (float) Math.pow(10, -45));
        System.out.println((float) Math.pow(10, 38));
        // Test string contains Float.MIN_VALUE, .MAX_VALUE, and +1 and -1 of them
        final String large = "affine 2 1.0 1.4E-45 3.4028235E38 " + (1.4E-45 + 1) + " " + (3.4028235E38 - 1) + " 1.0 2330.0 724.0 2330.0 724.0 1.0 2516.0 727.0 2516.0 727.0 1.0 2825.0 739.0 2825.0 739.0 1.0 3002.0002 736.0 3002.0 736.0 1.0 3281.001 736.0 3281.0 736.0 1.0 3482.001 742.0 3482.0 742.0 1.0 3725.001 745.0 3725.0 745.0 1.0 3983.0017 745.0 3983.0 745.0 1.0 4157.0024 739.0 4157.0 739.0 1.0 4325.002 751.0 4325.0 751.0 1.0 4390.9854 993.99994 4391.0 994.0 1.0 4219.9834 1020.9999 4220.0 1021.0 1.0 4018.983 1023.9999 4019.0 1024.0 1.0 3796.983 1017.9999 3797.0 1018.0 1.0 3586.9824 1020.9999 3587.0 1021.0 1.0 3289.981 1035.9998 3290.0 1036.0 1.0 3148.9802 1047.9998 3149.0 1048.0 1.0 2806.9795 1044.9998 2807.0 1045.0 1.0 2578.979 1047.9998 2579.0 1048.0 1.0 2392.9785 1047.9998 2393.0 1048.0 1.0 2371.969 1188.9998 2372.0 1189.0 1.0 2638.9685 1206.9998 2639.0 1207.0 1.0 2851.9688 1206.9998 2852.0 1207.0 1.0 3082.9697 1206.9998 3083.0 1207.0 1.0 3310.97 1206.9998 3311.0 1207.0 1.0 3496.9702 1206.9999 3497.0 1207.0 1.0 3691.9712 1198.0 3692.0 1198.0 1.0 3943.9717 1195.0 3944.0 1195.0 1.0 4123.9717 1195.0 4124.0 1195.0 1.0 4354.9727 1183.0 4355.0 1183.0 1.0 4423.965 1306.0 4424.0 1306.0 1.0 4231.9575 1410.9999 4232.0 1411.0 1.0 3850.957 1413.9999 3845.0 1414.0 1.0 3552.9067 1416.9998 3548.0 1417.0 1.0 3340.6917 1419.9998 3338.0 1420.0 1.0 3172.0889 1422.9995 3170.0 1423.0 1.0 3060.851 1422.9995 3059.0 1423.0 1.0 2823.295 1422.9995 2822.0 1423.0 1.0 2709.1892 1428.9995 2708.0 1429.0 1.0 2553.0269 1428.9995 2552.0 1429.0 1.0 2403.3828 1578.9995 2402.0 1579.0 1.0 2619.5566 1578.9995 2618.0 1579.0 1.0 2895.9922 1587.9996 2894.0 1588.0 1.0 3157.6333 1590.9996 3155.0 1591.0 1.0 3425.6125 1587.9998 3422.0 1588.0 1.0 3729.5825 1588.0 3725.0 1588.0 1.0 3894.6206 1588.0 3890.0 1588.0 1.0 4126.9688 1579.0 4124.0 1579.0 1.0 4348.025 1576.0 4346.0 1576.0 1.0 4452.919 1579.0002 4451.0 1579.0 1.0 4504.8438 1768.0002 4502.0 1768.0 1.0 4295.339 1798.0002 4292.0 1798.0 1.0 4091.9287 1807.0002 4088.0 1807.0 1.0 3792.756 1807.0001 3788.0 1807.0 1.0 3555.6267 1803.9998 3551.0 1804.0 1.0 3306.2095 1809.9998 3302.0 1810.0 1.0 3074.7664 1821.9998 3071.0 1822.0 1.0 2858.3079 1818.9995 2852.0 1819.0 1.0 2621.9727 1818.9995 2618.0 1819.0 1.0 2339.1707 1818.9995 2336.0 1819.0 1.0 2376.23 1998.9995 2372.0 1999.0 1.0 2643.9084 2001.9998 2639.0 2002.0 1.0 2926.463 2010.9998 2924.0 2011.0 1.0 3228.7146 2013.9998 3227.0 2014.0 1.0 3383.458 2016.9999 3380.0 2017.0 1.0 3621.5615 2011.0 3617.0 2011.0 1.0 3828.855 2011.0001 3824.0 2011.0 1.0 4059.7258 2008.0005 4055.0 2008.0 1.0 4179.616 2005.0005 4175.0 2005.0 1.0 4482.1504 1990.0006 4478.0 1990.0 1.0 4501.206 2215.0005 4496.0 2215.0 1.0 4339.3223 2221.0005 4334.0 2221.0 1.0 4009.5066 2224.0005 4004.0 2224.0 1.0 3718.4097 2233.0002 3707.0 2233.0 1.0 3385.523 2227.0002 3380.0 2227.0 1.0 3084.8071 2236.0 3080.0 2236.0 1.0 2850.8462 2235.9998 2843.0 2236.0 1.0 2618.0845 2229.9998 2606.0 2230.0 1.0 2520.2085 2232.9995 2510.0 2233.0 1.0 2305.331 2232.9995 2297.0 2233.0 1.0 2312.5972 2376.9995 2303.0 2380.0 1.0 2613.6096 2379.6626 2603.0 2380.0 1.0 2839.9526 2379.8416 2834.0 2380.0 1.0 2999.8557 2379.891 2993.0 2380.0 1.0 3206.6636 2379.927 3200.0 2380.0 1.0 3441.0208 2382.9553 3434.0 2383.0 1.0 3673.747 2382.9705 3665.0 2383.0 1.0 3853.3145 2382.979 3845.0 2383.0 1.0 4190.8413 2385.987 4184.0 2386.0 1.0 4487.301 2385.996 4481.0 2386.0 1.0 2296.2778 604.01086 2297.0 604.0 1.0 2536.1816 595.017 2537.0 595.0 1.0 2839.126 595.01697 2840.0 595.0 1.0 3118.0815 598.0175 3119.0 598.0 1.0 3364.12 601.01544 3365.0 601.0 1.0 3508.1118 601.01587 3509.0 601.0 1.0 3628.0789 598.01654 3638.0 598.0 1.0 3817.8489 598.016 3824.0 598.0 1.0 4019.6694 604.01337 4022.0 604.0 1.0 4278.367 583.0154 4280.0 583.0 1.0";
        final MovingLeastSquaresTransform2 m = new MovingLeastSquaresTransform2();
        // 
        Field fp = mpicbg.models.MovingLeastSquaresTransform2.class.getDeclaredField("p");
        fp.setAccessible(true);
        Field fq = mpicbg.models.MovingLeastSquaresTransform2.class.getDeclaredField("q");
        fq.setAccessible(true);
        Field fw = mpicbg.models.MovingLeastSquaresTransform2.class.getDeclaredField("w");
        fw.setAccessible(true);
        // Reference: Saalfeld's
        m.init(large);
        float a1 = m.getAlpha();
        float[][] p1 = (float[][]) fp.get(m);
        float[][] q1 = (float[][]) fq.get(m);
        float[] w1 = (float[]) fw.get(m);
        // To compare with: mine
        m.init2(large);
        float a2 = m.getAlpha();
        float[][] p2 = (float[][]) fp.get(m);
        float[][] q2 = (float[][]) fq.get(m);
        float[] w2 = (float[]) fw.get(m);
        // Test correctness: compare equality of all fields
        int nErrors = 0;
        if (a1 != a2) {
            System.out.println("ERROR with alpha: " + a1 + " != " + a2);
            ++nErrors;
        }
        for (int k = p1[0].length - 1; k > -1; --k) {
            for (int d = 0; d < p1.length; ++d) {
                if (p1[d][k] != p2[d][k]) {
                    System.out.println("ERROR with p: " + p1[d][k] + " != " + p2[d][k]);
                    ++nErrors;
                }
            }
            for (int d = 0; d < q1.length; ++d) {
                if (q1[d][k] != q2[d][k]) {
                    System.out.println("ERROR with q: " + q1[d][k] + " != " + q2[d][k]);
                    ++nErrors;
                }
            }
            if (w1[k] != w2[k]) {
                System.out.println("ERROR with w: " + w1[k] + " != " + w2[k]);
                ++nErrors;
            }
        }
        System.out.println("Number of differences between init1 and init2: " + nErrors);
        // Compare performance
        for (int k = 0; k < 10; ++k) {
            long t0 = System.currentTimeMillis();
            new MovingLeastSquaresTransform2().init(large);
            long t1 = System.currentTimeMillis();
            System.out.println("init1: " + (t1 - t0) + " ms");
        }
        for (int k = 0; k < 10; ++k) {
            long t0 = System.currentTimeMillis();
            new MovingLeastSquaresTransform2().init2(large);
            long t1 = System.currentTimeMillis();
            System.out.println("init2: " + (t1 - t0) + " ms");
        }
    } catch (Exception e) {
        e.printStackTrace();
    }
}
Also used : Field(java.lang.reflect.Field) MovingLeastSquaresTransform2(mpicbg.trakem2.transform.MovingLeastSquaresTransform2)

Example 5 with MovingLeastSquaresTransform2

use of mpicbg.trakem2.transform.MovingLeastSquaresTransform2 in project TrakEM2 by trakem2.

the class AlignTask method alignMultiLayerMosaicTask.

@SuppressWarnings({ "unchecked", "rawtypes" })
public static final void alignMultiLayerMosaicTask(final List<Layer> layerRange, final Patch nail, final Align.Param cp, final Align.ParamOptimize p, final Align.ParamOptimize pcp, final boolean tilesAreInPlaceIn, final boolean largestGraphOnlyIn, final boolean hideDisconnectedTilesIn, final boolean deleteDisconnectedTilesIn, final boolean deformIn) {
    /* register */
    final List<AbstractAffineTile2D<?>> allTiles = new ArrayList<AbstractAffineTile2D<?>>();
    final List<AbstractAffineTile2D<?>> allFixedTiles = new ArrayList<AbstractAffineTile2D<?>>();
    final List<AbstractAffineTile2D<?>> previousLayerTiles = new ArrayList<AbstractAffineTile2D<?>>();
    final HashMap<Patch, PointMatch> tileCenterPoints = new HashMap<Patch, PointMatch>();
    final Collection<Patch> fixedPatches = new HashSet<Patch>();
    if (null != nail)
        fixedPatches.add(nail);
    Layer previousLayer = null;
    for (final Layer layer : layerRange) {
        /* align all tiles in the layer */
        final List<Patch> patches = new ArrayList<Patch>();
        for (// ignore hidden tiles
        final Displayable a : // ignore hidden tiles
        layer.getDisplayables(Patch.class, true)) if (a instanceof Patch)
            patches.add((Patch) a);
        final List<AbstractAffineTile2D<?>> currentLayerTiles = new ArrayList<AbstractAffineTile2D<?>>();
        final List<AbstractAffineTile2D<?>> fixedTiles = new ArrayList<AbstractAffineTile2D<?>>();
        Align.tilesFromPatches(p, patches, fixedPatches, currentLayerTiles, fixedTiles);
        // Will consider graphs and hide/delete tiles when all cross-layer graphs are found.
        alignTiles(p, currentLayerTiles, fixedTiles, tilesAreInPlaceIn, false, false, false);
        if (Thread.currentThread().isInterrupted())
            return;
        /* connect to the previous layer */
        /* generate tiles with the cross-section model from the current layer tiles */
        /* ------------------------------------------------------------------------ */
        /* TODO step back and make tiles bare containers for a patch and a model such that by changing the model the tile can be reused */
        final HashMap<Patch, AbstractAffineTile2D<?>> currentLayerPatchTiles = new HashMap<Patch, AbstractAffineTile2D<?>>();
        for (final AbstractAffineTile2D<?> t : currentLayerTiles) currentLayerPatchTiles.put(t.getPatch(), t);
        final List<AbstractAffineTile2D<?>> csCurrentLayerTiles = new ArrayList<AbstractAffineTile2D<?>>();
        final Set<AbstractAffineTile2D<?>> csFixedTiles = new HashSet<AbstractAffineTile2D<?>>();
        Align.tilesFromPatches(cp, patches, fixedPatches, csCurrentLayerTiles, csFixedTiles);
        final HashMap<Tile<?>, AbstractAffineTile2D<?>> tileTiles = new HashMap<Tile<?>, AbstractAffineTile2D<?>>();
        for (final AbstractAffineTile2D<?> t : csCurrentLayerTiles) tileTiles.put(currentLayerPatchTiles.get(t.getPatch()), t);
        for (final AbstractAffineTile2D<?> t : currentLayerTiles) {
            final AbstractAffineTile2D<?> csLayerTile = tileTiles.get(t);
            csLayerTile.addMatches(t.getMatches());
            for (final Tile<?> ct : t.getConnectedTiles()) csLayerTile.addConnectedTile(tileTiles.get(ct));
        }
        /* add a fixed tile only if there was a Patch selected */
        allFixedTiles.addAll(csFixedTiles);
        /* first, align connected graphs to each other */
        /* graphs in the current layer */
        final List<Set<Tile<?>>> currentLayerGraphs = AbstractAffineTile2D.identifyConnectedGraphs(csCurrentLayerTiles);
        if (Thread.currentThread().isInterrupted())
            return;
        // /* TODO just for visualization */
        // for ( final Set< Tile< ? > > graph : currentLayerGraphs )
        // {
        // Display.getFront().getSelection().clear();
        // Display.getFront().setLayer( ( ( AbstractAffineTile2D< ? > )graph.iterator().next() ).getPatch().getLayer() );
        // 
        // for ( final Tile< ? > tile : graph )
        // {
        // Display.getFront().getSelection().add( ( ( AbstractAffineTile2D< ? > )tile ).getPatch() );
        // Display.repaint();
        // }
        // Utils.showMessage( "OK" );
        // }
        /* graphs from the whole system that are present in the previous layer */
        final List<Set<Tile<?>>> graphs = AbstractAffineTile2D.identifyConnectedGraphs(allTiles);
        final HashMap<Set<Tile<?>>, Set<Tile<?>>> graphGraphs = new HashMap<Set<Tile<?>>, Set<Tile<?>>>();
        for (final Set<Tile<?>> graph : graphs) {
            if (Thread.currentThread().isInterrupted())
                return;
            final Set<Tile<?>> previousLayerGraph = new HashSet<Tile<?>>();
            for (final Tile<?> tile : previousLayerTiles) {
                if (graph.contains(tile)) {
                    graphGraphs.put(graph, previousLayerGraph);
                    previousLayerGraph.add(tile);
                }
            }
        }
        final Collection<Set<Tile<?>>> previousLayerGraphs = graphGraphs.values();
        // /* TODO just for visualization */
        // for ( final Set< Tile< ? > > graph : previousLayerGraphs )
        // {
        // Display.getFront().getSelection().clear();
        // Display.getFront().setLayer( ( ( AbstractAffineTile2D< ? > )graph.iterator().next() ).getPatch().getLayer() );
        // 
        // for ( final Tile< ? > tile : graph )
        // {
        // Display.getFront().getSelection().add( ( ( AbstractAffineTile2D< ? > )tile ).getPatch() );
        // Display.repaint();
        // }
        // Utils.showMessage( "OK" );
        // }
        /* generate snapshots of the graphs and preregister them using the parameters defined in cp */
        final List<AbstractAffineTile2D<?>[]> crossLayerTilePairs = new ArrayList<AbstractAffineTile2D<?>[]>();
        for (final Set<Tile<?>> currentLayerGraph : currentLayerGraphs) {
            for (final Set<Tile<?>> previousLayerGraph : previousLayerGraphs) {
                if (Thread.currentThread().isInterrupted())
                    return;
                alignGraphs(cp, layer, previousLayer, currentLayerGraph, previousLayerGraph);
                /* TODO this is pointless data shuffling just for type incompatibility---fix this at the root */
                final ArrayList<AbstractAffineTile2D<?>> previousLayerGraphTiles = new ArrayList<AbstractAffineTile2D<?>>();
                previousLayerGraphTiles.addAll((Set) previousLayerGraph);
                final ArrayList<AbstractAffineTile2D<?>> currentLayerGraphTiles = new ArrayList<AbstractAffineTile2D<?>>();
                currentLayerGraphTiles.addAll((Set) currentLayerGraph);
                AbstractAffineTile2D.pairOverlappingTiles(previousLayerGraphTiles, currentLayerGraphTiles, crossLayerTilePairs);
            }
        }
        /* ------------------------------------------------------------------------ */
        /* this is without the affine/rigid approximation per graph */
        // AbstractAffineTile2D.pairTiles( previousLayerTiles, csCurrentLayerTiles, crossLayerTilePairs );
        Align.connectTilePairs(cp, csCurrentLayerTiles, crossLayerTilePairs, Runtime.getRuntime().availableProcessors());
        if (Thread.currentThread().isInterrupted())
            return;
        // for ( final AbstractAffineTile2D< ? >[] tilePair : crossLayerTilePairs )
        // {
        // Display.getFront().setLayer( tilePair[ 0 ].getPatch().getLayer() );
        // Display.getFront().getSelection().clear();
        // Display.getFront().getSelection().add( tilePair[ 0 ].getPatch() );
        // Display.getFront().getSelection().add( tilePair[ 1 ].getPatch() );
        // 
        // Utils.showMessage( "1: OK?" );
        // 
        // Display.getFront().setLayer( tilePair[ 1 ].getPatch().getLayer() );
        // Display.getFront().getSelection().clear();
        // Display.getFront().getSelection().add( tilePair[ 0 ].getPatch() );
        // Display.getFront().getSelection().add( tilePair[ 1 ].getPatch() );
        // 
        // Utils.showMessage( "2: OK?" );
        // }
        /* prepare the next loop */
        allTiles.addAll(csCurrentLayerTiles);
        previousLayerTiles.clear();
        previousLayerTiles.addAll(csCurrentLayerTiles);
        /* optimize */
        Align.optimizeTileConfiguration(pcp, allTiles, allFixedTiles);
        if (Thread.currentThread().isInterrupted())
            return;
        for (final AbstractAffineTile2D<?> t : allTiles) t.getPatch().setAffineTransform(t.getModel().createAffine());
        previousLayer = layer;
    }
    final List<Set<Tile<?>>> graphs = AbstractAffineTile2D.identifyConnectedGraphs(allTiles);
    final List<AbstractAffineTile2D<?>> interestingTiles = new ArrayList<AbstractAffineTile2D<?>>();
    if (largestGraphOnlyIn && (hideDisconnectedTilesIn || deleteDisconnectedTilesIn)) {
        if (Thread.currentThread().isInterrupted())
            return;
        /* find largest graph. */
        Set<Tile<?>> largestGraph = null;
        for (final Set<Tile<?>> graph : graphs) if (largestGraph == null || largestGraph.size() < graph.size())
            largestGraph = graph;
        final Set<AbstractAffineTile2D<?>> tiles_to_keep = new HashSet<AbstractAffineTile2D<?>>();
        for (final Tile<?> t : largestGraph) tiles_to_keep.add((AbstractAffineTile2D<?>) t);
        if (hideDisconnectedTilesIn)
            for (final AbstractAffineTile2D<?> t : allTiles) if (!tiles_to_keep.contains(t))
                t.getPatch().setVisible(false);
        if (deleteDisconnectedTilesIn)
            for (final AbstractAffineTile2D<?> t : allTiles) if (!tiles_to_keep.contains(t))
                t.getPatch().remove(false);
        interestingTiles.addAll(tiles_to_keep);
    } else
        interestingTiles.addAll(allTiles);
    if (deformIn) {
        /* ############################################ */
        /* experimental: use the center points of all tiles to define a MLS deformation from the pure intra-layer registration to the globally optimal */
        Utils.log("deforming...");
        /* store the center location of each single tile for later deformation */
        for (final AbstractAffineTile2D<?> t : interestingTiles) {
            final double[] c = new double[] { t.getWidth() / 2.0, t.getHeight() / 2.0 };
            t.getModel().applyInPlace(c);
            final Point q = new Point(c);
            tileCenterPoints.put(t.getPatch(), new PointMatch(q.clone(), q));
        }
        for (final Layer layer : layerRange) {
            Utils.log("layer" + layer);
            if (Thread.currentThread().isInterrupted())
                return;
            /* again, align all tiles in the layer */
            final List<Patch> patches = new ArrayList<Patch>();
            for (final Displayable a : layer.getDisplayables(Patch.class)) if (a instanceof Patch)
                patches.add((Patch) a);
            final List<AbstractAffineTile2D<?>> currentLayerTiles = new ArrayList<AbstractAffineTile2D<?>>();
            final List<AbstractAffineTile2D<?>> fixedTiles = new ArrayList<AbstractAffineTile2D<?>>();
            Align.tilesFromPatches(p, patches, fixedPatches, currentLayerTiles, fixedTiles);
            /* add a fixed tile only if there was a Patch selected */
            allFixedTiles.addAll(fixedTiles);
            // will consider graphs and hide/delete tiles when all cross-layer graphs are found
            alignTiles(p, currentLayerTiles, fixedTiles, true, false, false, false);
            /* for each independent graph do an independent transform */
            final List<Set<Tile<?>>> currentLayerGraphs = AbstractAffineTile2D.identifyConnectedGraphs(currentLayerTiles);
            for (final Set<Tile<?>> graph : currentLayerGraphs) {
                /* update the tile-center pointmatches */
                final Collection<PointMatch> matches = new ArrayList<PointMatch>();
                final Collection<AbstractAffineTile2D<?>> toBeDeformedTiles = new ArrayList<AbstractAffineTile2D<?>>();
                for (final AbstractAffineTile2D<?> t : (Collection<AbstractAffineTile2D<?>>) (Collection) graph) {
                    final PointMatch pm = tileCenterPoints.get(t.getPatch());
                    if (pm == null)
                        continue;
                    final double[] pl = pm.getP1().getL();
                    pl[0] = t.getWidth() / 2.0;
                    pl[1] = t.getHeight() / 2.0;
                    t.getModel().applyInPlace(pl);
                    matches.add(pm);
                    toBeDeformedTiles.add(t);
                }
                for (final AbstractAffineTile2D<?> t : toBeDeformedTiles) {
                    if (Thread.currentThread().isInterrupted())
                        return;
                    try {
                        final Patch patch = t.getPatch();
                        final Rectangle pbox = patch.getCoordinateTransformBoundingBox();
                        final AffineTransform pat = new AffineTransform();
                        pat.translate(-pbox.x, -pbox.y);
                        pat.preConcatenate(patch.getAffineTransform());
                        final mpicbg.trakem2.transform.AffineModel2D toWorld = new mpicbg.trakem2.transform.AffineModel2D();
                        toWorld.set(pat);
                        final MovingLeastSquaresTransform2 mlst = Align.createMLST(matches, 1.0f);
                        final CoordinateTransformList<CoordinateTransform> ctl = new CoordinateTransformList<CoordinateTransform>();
                        ctl.add(toWorld);
                        ctl.add(mlst);
                        ctl.add(toWorld.createInverse());
                        patch.appendCoordinateTransform(ctl);
                        patch.getProject().getLoader().regenerateMipMaps(patch);
                    } catch (final Exception e) {
                        e.printStackTrace();
                    }
                }
            }
        }
    }
    layerRange.get(0).getParent().setMinimumDimensions();
    IJ.log("Done: register multi-layer mosaic.");
    return;
}
Also used : Set(java.util.Set) HashSet(java.util.HashSet) LayerSet(ini.trakem2.display.LayerSet) HashMap(java.util.HashMap) CoordinateTransformList(mpicbg.trakem2.transform.CoordinateTransformList) ArrayList(java.util.ArrayList) Rectangle(java.awt.Rectangle) AbstractAffineModel2D(mpicbg.models.AbstractAffineModel2D) AffineModel2D(mpicbg.models.AffineModel2D) HashSet(java.util.HashSet) Displayable(ini.trakem2.display.Displayable) Tile(mpicbg.models.Tile) Point(mpicbg.models.Point) Layer(ini.trakem2.display.Layer) NotEnoughDataPointsException(mpicbg.models.NotEnoughDataPointsException) NoninvertibleModelException(mpicbg.models.NoninvertibleModelException) NoninvertibleTransformException(java.awt.geom.NoninvertibleTransformException) PointMatch(mpicbg.models.PointMatch) MovingLeastSquaresTransform2(mpicbg.trakem2.transform.MovingLeastSquaresTransform2) Collection(java.util.Collection) AffineTransform(java.awt.geom.AffineTransform) Patch(ini.trakem2.display.Patch) CoordinateTransform(mpicbg.trakem2.transform.CoordinateTransform) InvertibleCoordinateTransform(mpicbg.trakem2.transform.InvertibleCoordinateTransform)

Aggregations

MovingLeastSquaresTransform2 (mpicbg.trakem2.transform.MovingLeastSquaresTransform2)5 ArrayList (java.util.ArrayList)3 AbstractAffineModel2D (mpicbg.models.AbstractAffineModel2D)3 AffineModel2D (mpicbg.models.AffineModel2D)3 Point (mpicbg.models.Point)3 PointMatch (mpicbg.models.PointMatch)3 Layer (ini.trakem2.display.Layer)2 LayerSet (ini.trakem2.display.LayerSet)2 Patch (ini.trakem2.display.Patch)2 NotEnoughDataPointsException (mpicbg.models.NotEnoughDataPointsException)2 Tile (mpicbg.models.Tile)2 CoordinateTransform (mpicbg.trakem2.transform.CoordinateTransform)2 Displayable (ini.trakem2.display.Displayable)1 VectorData (ini.trakem2.display.VectorData)1 Rectangle (java.awt.Rectangle)1 AffineTransform (java.awt.geom.AffineTransform)1 Area (java.awt.geom.Area)1 NoninvertibleTransformException (java.awt.geom.NoninvertibleTransformException)1 Field (java.lang.reflect.Field)1 Collection (java.util.Collection)1