Examples with Filter ini.trakem2.utils.Filter used on opensource projects

Example 6 with Filter

use of ini.trakem2.utils.Filter 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.");
}

Example 7 with Filter

use of ini.trakem2.utils.Filter in project TrakEM2 by trakem2.

the class ElasticMontage method exec.

@SuppressWarnings("deprecation")
public final void exec(final Param param, final List<Patch> patches, final Set<Patch> fixedPatches) throws Exception {
    /* free memory */
    patches.get(0).getProject().getLoader().releaseAll();
    /* create tiles and models for all patches */
    final ArrayList<AbstractAffineTile2D<?>> tiles = new ArrayList<AbstractAffineTile2D<?>>();
    final ArrayList<AbstractAffineTile2D<?>> fixedTiles = new ArrayList<AbstractAffineTile2D<?>>();
    Align.tilesFromPatches(param.po, patches, fixedPatches, tiles, fixedTiles);
    if (!param.isAligned) {
        Align.alignTiles(param.po, tiles, fixedTiles, param.tilesAreInPlace, param.maxNumThreads);
        /* Apply the estimated affine transform to patches */
        for (final AbstractAffineTile2D<?> t : tiles) t.getPatch().setAffineTransform(t.createAffine());
        Display.update();
    }
    /* generate tile pairs for all by now overlapping tiles */
    final ArrayList<AbstractAffineTile2D<?>[]> tilePairs = new ArrayList<AbstractAffineTile2D<?>[]>();
    AbstractAffineTile2D.pairOverlappingTiles(tiles, tilePairs);
    /* check if there was any pair */
    if (tilePairs.size() == 0) {
        Utils.log("Elastic montage could not find any overlapping patches after pre-montaging.");
        return;
    }
    Utils.log(tilePairs.size() + " pairs of patches will be block-matched...");
    /* make pairwise global models local */
    final ArrayList<Triple<AbstractAffineTile2D<?>, AbstractAffineTile2D<?>, InvertibleCoordinateTransform>> pairs = new ArrayList<Triple<AbstractAffineTile2D<?>, AbstractAffineTile2D<?>, InvertibleCoordinateTransform>>();
    /*
		 * The following casting madness is necessary to get this code compiled
		 * with Sun/Oracle Java 6 which otherwise generates an inconvertible
		 * type exception.
		 *
		 * TODO Remove as soon as this bug is fixed in Sun/Oracle javac.
		 */
    for (final AbstractAffineTile2D<?>[] pair : tilePairs) {
        final AbstractAffineModel2D<?> m;
        switch(param.po.desiredModelIndex) {
            case 0:
                final TranslationModel2D t = (TranslationModel2D) (Object) pair[1].getModel().createInverse();
                t.concatenate((TranslationModel2D) (Object) pair[0].getModel());
                m = t;
                break;
            case 1:
                final RigidModel2D r = (RigidModel2D) (Object) pair[1].getModel().createInverse();
                r.concatenate((RigidModel2D) (Object) pair[0].getModel());
                m = r;
                break;
            case 2:
                final SimilarityModel2D s = (SimilarityModel2D) (Object) pair[1].getModel().createInverse();
                s.concatenate((SimilarityModel2D) (Object) pair[0].getModel());
                m = s;
                break;
            case 3:
                final AffineModel2D a = (AffineModel2D) (Object) pair[1].getModel().createInverse();
                a.concatenate((AffineModel2D) (Object) pair[0].getModel());
                m = a;
                break;
            default:
                m = null;
        }
        pairs.add(new Triple<AbstractAffineTile2D<?>, AbstractAffineTile2D<?>, InvertibleCoordinateTransform>(pair[0], pair[1], m));
    }
    /* Elastic alignment */
    /* Initialization */
    final double springTriangleHeightTwice = 2 * Math.sqrt(0.75 * param.springLengthSpringMesh * param.springLengthSpringMesh);
    final ArrayList<SpringMesh> meshes = new ArrayList<SpringMesh>(tiles.size());
    final HashMap<AbstractAffineTile2D<?>, SpringMesh> tileMeshMap = new HashMap<AbstractAffineTile2D<?>, SpringMesh>();
    for (final AbstractAffineTile2D<?> tile : tiles) {
        final double w = tile.getWidth();
        final double h = tile.getHeight();
        final int numX = Math.max(2, (int) Math.ceil(w / param.springLengthSpringMesh) + 1);
        final int numY = Math.max(2, (int) Math.ceil(h / springTriangleHeightTwice) + 1);
        final double wMesh = (numX - 1) * param.springLengthSpringMesh;
        final double hMesh = (numY - 1) * springTriangleHeightTwice;
        final SpringMesh mesh = new SpringMesh(numX, numY, wMesh, hMesh, param.stiffnessSpringMesh, param.maxStretchSpringMesh * param.bmScale, param.dampSpringMesh);
        meshes.add(mesh);
        tileMeshMap.put(tile, mesh);
    }
    // final int blockRadius = Math.max( 32, Util.roundPos( param.springLengthSpringMesh / 2 ) );
    final int blockRadius = Math.max(Util.roundPos(16 / param.bmScale), param.bmBlockRadius);
    /**
     * TODO set this something more than the largest error by the approximate model
     */
    final int searchRadius = param.bmSearchRadius;
    final AbstractModel<?> localSmoothnessFilterModel = mpicbg.trakem2.align.Util.createModel(param.bmLocalModelIndex);
    for (final Triple<AbstractAffineTile2D<?>, AbstractAffineTile2D<?>, InvertibleCoordinateTransform> pair : pairs) {
        final AbstractAffineTile2D<?> t1 = pair.a;
        final AbstractAffineTile2D<?> t2 = pair.b;
        final SpringMesh m1 = tileMeshMap.get(t1);
        final SpringMesh m2 = tileMeshMap.get(t2);
        final ArrayList<PointMatch> pm12 = new ArrayList<PointMatch>();
        final ArrayList<PointMatch> pm21 = new ArrayList<PointMatch>();
        final ArrayList<Vertex> v1 = m1.getVertices();
        final ArrayList<Vertex> v2 = m2.getVertices();
        final String patchName1 = patchName(t1.getPatch());
        final String patchName2 = patchName(t2.getPatch());
        final PatchImage pi1 = t1.getPatch().createTransformedImage();
        if (pi1 == null) {
            Utils.log("Patch `" + patchName1 + "' failed generating a transformed image.  Skipping...");
            continue;
        }
        final PatchImage pi2 = t2.getPatch().createTransformedImage();
        if (pi2 == null) {
            Utils.log("Patch `" + patchName2 + "' failed generating a transformed image.  Skipping...");
            continue;
        }
        final FloatProcessor fp1 = (FloatProcessor) pi1.target.convertToFloat();
        final ByteProcessor mask1 = pi1.getMask();
        final FloatProcessor fpMask1 = mask1 == null ? null : scaleByte(mask1);
        final FloatProcessor fp2 = (FloatProcessor) pi2.target.convertToFloat();
        final ByteProcessor mask2 = pi2.getMask();
        final FloatProcessor fpMask2 = mask2 == null ? null : scaleByte(mask2);
        if (!fixedTiles.contains(t1)) {
            BlockMatching.matchByMaximalPMCC(fp1, fp2, fpMask1, fpMask2, param.bmScale, pair.c, blockRadius, blockRadius, searchRadius, searchRadius, param.bmMinR, param.bmRodR, param.bmMaxCurvatureR, v1, pm12, new ErrorStatistic(1));
            if (param.bmUseLocalSmoothnessFilter) {
                Utils.log("`" + patchName1 + "' > `" + patchName2 + "': found " + pm12.size() + " correspondence candidates.");
                localSmoothnessFilterModel.localSmoothnessFilter(pm12, pm12, param.bmLocalRegionSigma, param.bmMaxLocalEpsilon, param.bmMaxLocalTrust);
                Utils.log("`" + patchName1 + "' > `" + patchName2 + "': " + pm12.size() + " candidates passed local smoothness filter.");
            } else {
                Utils.log("`" + patchName1 + "' > `" + patchName2 + "': found " + pm12.size() + " correspondences.");
            }
        } else {
            Utils.log("Skipping fixed patch `" + patchName1 + "'.");
        }
        if (!fixedTiles.contains(t2)) {
            BlockMatching.matchByMaximalPMCC(fp2, fp1, fpMask2, fpMask1, param.bmScale, pair.c.createInverse(), blockRadius, blockRadius, searchRadius, searchRadius, param.bmMinR, param.bmRodR, param.bmMaxCurvatureR, v2, pm21, new ErrorStatistic(1));
            if (param.bmUseLocalSmoothnessFilter) {
                Utils.log("`" + patchName1 + "' < `" + patchName2 + "': found " + pm21.size() + " correspondence candidates.");
                localSmoothnessFilterModel.localSmoothnessFilter(pm21, pm21, param.bmLocalRegionSigma, param.bmMaxLocalEpsilon, param.bmMaxLocalTrust);
                Utils.log("`" + patchName1 + "' < `" + patchName2 + "': " + pm21.size() + " candidates passed local smoothness filter.");
            } else {
                Utils.log("`" + patchName1 + "' < `" + patchName2 + "': found " + pm21.size() + " correspondences.");
            }
        } else {
            Utils.log("Skipping fixed patch `" + patchName2 + "'.");
        }
        for (final PointMatch pm : pm12) {
            final Vertex p1 = (Vertex) pm.getP1();
            final Vertex p2 = new Vertex(pm.getP2());
            p1.addSpring(p2, new Spring(0, 1.0f));
            m2.addPassiveVertex(p2);
        }
        for (final PointMatch pm : pm21) {
            final Vertex p1 = (Vertex) pm.getP1();
            final Vertex p2 = new Vertex(pm.getP2());
            p1.addSpring(p2, new Spring(0, 1.0f));
            m1.addPassiveVertex(p2);
        }
    }
    /* initialize */
    for (final Map.Entry<AbstractAffineTile2D<?>, SpringMesh> entry : tileMeshMap.entrySet()) entry.getValue().init(entry.getKey().getModel());
    /* optimize the meshes */
    try {
        final long t0 = System.currentTimeMillis();
        IJ.log("Optimizing spring meshes...");
        if (param.useLegacyOptimizer) {
            Utils.log("  ...using legacy optimizer...");
            SpringMesh.optimizeMeshes2(meshes, param.po.maxEpsilon, param.maxIterationsSpringMesh, param.maxPlateauwidthSpringMesh, param.visualize);
        } else {
            SpringMesh.optimizeMeshes(meshes, param.po.maxEpsilon, param.maxIterationsSpringMesh, param.maxPlateauwidthSpringMesh, param.visualize);
        }
        IJ.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;
    }
    /* apply */
    for (final Map.Entry<AbstractAffineTile2D<?>, SpringMesh> entry : tileMeshMap.entrySet()) {
        final AbstractAffineTile2D<?> tile = entry.getKey();
        if (!fixedTiles.contains(tile)) {
            final Patch patch = tile.getPatch();
            final SpringMesh mesh = entry.getValue();
            final Set<PointMatch> matches = mesh.getVA().keySet();
            Rectangle box = patch.getCoordinateTransformBoundingBox();
            /* compensate for existing coordinate transform bounding box */
            for (final PointMatch pm : matches) {
                final Point p1 = pm.getP1();
                final double[] l = p1.getL();
                l[0] += box.x;
                l[1] += box.y;
            }
            final ThinPlateSplineTransform mlt = ElasticLayerAlignment.makeTPS(matches);
            patch.appendCoordinateTransform(mlt);
            box = patch.getCoordinateTransformBoundingBox();
            patch.getAffineTransform().setToTranslation(box.x, box.y);
            patch.updateInDatabase("transform");
            patch.updateBucket();
            patch.updateMipMaps();
        }
    }
    Utils.log("Done.");
}

Example 8 with Filter

use of ini.trakem2.utils.Filter in project TrakEM2 by trakem2.

the class AlignLayersTask method alignLayersNonLinearlyJob.

public static final void alignLayersNonLinearlyJob(final LayerSet layerSet, final int first, final int last, final boolean propagateTransform, final Rectangle fov, final Filter<Patch> filter) {
    // will reverse order if necessary
    final List<Layer> layerRange = layerSet.getLayers(first, last);
    final Align.Param p = Align.param.clone();
    // Remove all empty layers
    for (final Iterator<Layer> it = layerRange.iterator(); it.hasNext(); ) {
        if (!it.next().contains(Patch.class, true)) {
            it.remove();
        }
    }
    if (0 == layerRange.size()) {
        Utils.log("No layers in range show any images!");
        return;
    }
    /* do not work if there is only one layer selected */
    if (layerRange.size() < 2)
        return;
    final List<Patch> all = new ArrayList<Patch>();
    for (final Layer la : layerRange) {
        for (final Patch patch : la.getAll(Patch.class)) {
            if (null != filter && !filter.accept(patch))
                continue;
            all.add(patch);
        }
    }
    AlignTask.transformPatchesAndVectorData(all, new Runnable() {

        @Override
        public void run() {
            // ///
            final Loader loader = layerSet.getProject().getLoader();
            // Not concurrent safe! So two copies, one per layer and Thread:
            final SIFT ijSIFT1 = new SIFT(new FloatArray2DSIFT(p.sift));
            final SIFT ijSIFT2 = new SIFT(new FloatArray2DSIFT(p.sift));
            final Collection<Feature> features1 = new ArrayList<Feature>();
            final Collection<Feature> features2 = new ArrayList<Feature>();
            final List<PointMatch> candidates = new ArrayList<PointMatch>();
            final List<PointMatch> inliers = new ArrayList<PointMatch>();
            final int n_proc = Runtime.getRuntime().availableProcessors() > 1 ? 2 : 1;
            final ExecutorService exec = Utils.newFixedThreadPool(n_proc, "alignLayersNonLinearly");
            List<Patch> previousPatches = null;
            int s = 0;
            for (int i = 1; i < layerRange.size(); ++i) {
                if (Thread.currentThread().isInterrupted())
                    break;
                final Layer layer1 = layerRange.get(i - 1);
                final Layer layer2 = layerRange.get(i);
                final long t0 = System.currentTimeMillis();
                features1.clear();
                features2.clear();
                final Rectangle box1 = null == fov ? layer1.getMinimalBoundingBox(Patch.class, true) : fov;
                final Rectangle box2 = null == fov ? layer2.getMinimalBoundingBox(Patch.class, true) : fov;
                /* calculate the common scale factor for both flat images */
                final double scale = Math.min(1.0f, (double) p.sift.maxOctaveSize / (double) Math.max(box1.width, Math.max(box1.height, Math.max(box2.width, box2.height))));
                final List<Patch> patches1;
                if (null == previousPatches) {
                    patches1 = layer1.getAll(Patch.class);
                    if (null != filter) {
                        for (final Iterator<Patch> it = patches1.iterator(); it.hasNext(); ) {
                            if (!filter.accept(it.next()))
                                it.remove();
                        }
                    }
                } else {
                    patches1 = previousPatches;
                }
                final List<Patch> patches2 = layer2.getAll(Patch.class);
                if (null != filter) {
                    for (final Iterator<Patch> it = patches2.iterator(); it.hasNext(); ) {
                        if (!filter.accept(it.next()))
                            it.remove();
                    }
                }
                final Future<ImageProcessor> fu1 = exec.submit(new Callable<ImageProcessor>() {

                    @Override
                    public ImageProcessor call() {
                        final ImageProcessor ip1 = loader.getFlatImage(layer1, box1, scale, 0xffffffff, ImagePlus.GRAY8, Patch.class, patches1, true).getProcessor();
                        ijSIFT1.extractFeatures(ip1, features1);
                        Utils.log(features1.size() + " features extracted in layer \"" + layer1.getTitle() + "\" (took " + (System.currentTimeMillis() - t0) + " ms).");
                        return ip1;
                    }
                });
                final Future<ImageProcessor> fu2 = exec.submit(new Callable<ImageProcessor>() {

                    @Override
                    public ImageProcessor call() {
                        final ImageProcessor ip2 = loader.getFlatImage(layer2, box2, scale, 0xffffffff, ImagePlus.GRAY8, Patch.class, patches2, true).getProcessor();
                        ijSIFT2.extractFeatures(ip2, features2);
                        Utils.log(features2.size() + " features extracted in layer \"" + layer2.getTitle() + "\" (took " + (System.currentTimeMillis() - t0) + " ms).");
                        return ip2;
                    }
                });
                final ImageProcessor ip1, ip2;
                try {
                    ip1 = fu1.get();
                    ip2 = fu2.get();
                } catch (final Exception e) {
                    IJError.print(e);
                    return;
                }
                if (features1.size() > 0 && features2.size() > 0) {
                    final long t1 = System.currentTimeMillis();
                    candidates.clear();
                    FeatureTransform.matchFeatures(features2, features1, candidates, p.rod);
                    final AbstractAffineModel2D<?> model;
                    switch(p.expectedModelIndex) {
                        case 0:
                            model = new TranslationModel2D();
                            break;
                        case 1:
                            model = new RigidModel2D();
                            break;
                        case 2:
                            model = new SimilarityModel2D();
                            break;
                        case 3:
                            model = new AffineModel2D();
                            break;
                        default:
                            return;
                    }
                    boolean modelFound;
                    boolean again = false;
                    try {
                        do {
                            again = false;
                            modelFound = model.filterRansac(candidates, inliers, 1000, p.maxEpsilon, p.minInlierRatio, p.minNumInliers, 3);
                            if (modelFound && p.rejectIdentity) {
                                final ArrayList<Point> points = new ArrayList<Point>();
                                PointMatch.sourcePoints(inliers, points);
                                if (Transforms.isIdentity(model, points, p.identityTolerance)) {
                                    IJ.log("Identity transform for " + inliers.size() + " matches rejected.");
                                    candidates.removeAll(inliers);
                                    inliers.clear();
                                    again = true;
                                }
                            }
                        } while (again);
                    } catch (final NotEnoughDataPointsException e) {
                        modelFound = false;
                    }
                    if (modelFound) {
                        IJ.log("Model found for layer \"" + layer2.getTitle() + "\" and its predecessor:\n  correspondences  " + inliers.size() + " of " + candidates.size() + "\n  average residual error  " + (model.getCost() / scale) + " px\n  took " + (System.currentTimeMillis() - t1) + " ms");
                        final ImagePlus imp1 = new ImagePlus("target", ip1);
                        final ImagePlus imp2 = new ImagePlus("source", ip2);
                        final List<Point> sourcePoints = new ArrayList<Point>();
                        final List<Point> targetPoints = new ArrayList<Point>();
                        PointMatch.sourcePoints(inliers, sourcePoints);
                        PointMatch.targetPoints(inliers, targetPoints);
                        imp2.setRoi(Util.pointsToPointRoi(sourcePoints));
                        imp1.setRoi(Util.pointsToPointRoi(targetPoints));
                        final ImageProcessor mask1 = ip1.duplicate();
                        mask1.threshold(1);
                        final ImageProcessor mask2 = ip2.duplicate();
                        mask2.threshold(1);
                        final Transformation warp = bUnwarpJ_.computeTransformationBatch(imp2, imp1, mask2, mask1, elasticParam);
                        final CubicBSplineTransform transf = new CubicBSplineTransform();
                        transf.set(warp.getIntervals(), warp.getDirectDeformationCoefficientsX(), warp.getDirectDeformationCoefficientsY(), imp2.getWidth(), imp2.getHeight());
                        final ArrayList<Future<?>> fus = new ArrayList<Future<?>>();
                        // Transform desired patches only
                        for (final Patch patch : patches2) {
                            try {
                                final Rectangle pbox = patch.getCoordinateTransformBoundingBox();
                                final AffineTransform at = patch.getAffineTransform();
                                final AffineTransform pat = new AffineTransform();
                                pat.scale(scale, scale);
                                pat.translate(-box2.x, -box2.y);
                                pat.concatenate(at);
                                pat.translate(-pbox.x, -pbox.y);
                                final mpicbg.trakem2.transform.AffineModel2D toWorld = new mpicbg.trakem2.transform.AffineModel2D();
                                toWorld.set(pat);
                                final CoordinateTransformList<CoordinateTransform> ctl = new CoordinateTransformList<CoordinateTransform>();
                                // move the patch into the global space where bUnwarpJ calculated the transformation
                                ctl.add(toWorld);
                                // Apply non-linear transformation
                                ctl.add(transf);
                                // move it back
                                ctl.add(toWorld.createInverse());
                                patch.appendCoordinateTransform(ctl);
                                fus.add(patch.updateMipMaps());
                                // Compensate for offset between boxes
                                final AffineTransform offset = new AffineTransform();
                                offset.translate(box1.x - box2.x, box1.y - box2.y);
                                offset.concatenate(at);
                                patch.setAffineTransform(offset);
                            } catch (final Exception e) {
                                e.printStackTrace();
                            }
                        }
                        // await regeneration of all mipmaps
                        Utils.wait(fus);
                        Display.repaint(layer2);
                    } else
                        IJ.log("No model found for layer \"" + layer2.getTitle() + "\" and its predecessor:\n  correspondence candidates  " + candidates.size() + "\n  took " + (System.currentTimeMillis() - s) + " ms");
                }
                IJ.showProgress(++s, layerRange.size());
                // for next iteration
                previousPatches = patches2;
            }
            exec.shutdown();
            if (propagateTransform)
                Utils.log("Propagation not implemented yet for non-linear layer alignment.");
        /* // CANNOT be done (at least not trivially:
		 * //an appropriate "scale" cannot be computed, and the box2 is part of the spline computation.
		if ( propagateTransform && null != lastTransform )
		{
			for (final Layer la : l.getParent().getLayers(last > first ? last +1 : first -1, last > first ? l.getParent().size() -1 : 0)) {
				// Transform visible patches only
				final Rectangle box2 = la.getMinimalBoundingBox( Patch.class, true );
				for ( final Displayable disp : la.getDisplayables( Patch.class, true ) )
				{
					// ...
				}
			}
		}
		*/
        }
    });
// end of transformPatchesAndVectorData
}

Example 9 with Filter

use of ini.trakem2.utils.Filter in project TrakEM2 by trakem2.

the class AlignLayersTask method alignLayers.

public static final void alignLayers(final Layer l, final Rectangle fov) throws Exception {
    final List<Layer> layers = l.getParent().getLayers();
    final String[] layerTitles = new String[layers.size()];
    final String[] noneAndLayerTitles = new String[layers.size() + 1];
    noneAndLayerTitles[0] = "*None*";
    for (int i = 0; i < layers.size(); ++i) {
        layerTitles[i] = l.getProject().findLayerThing(layers.get(i)).toString();
        noneAndLayerTitles[i + 1] = layerTitles[i];
    }
    // Param p = Align.param;
    // Align.param.sift.maxOctaveSize = 1024;
    final GenericDialog gd = new GenericDialog("Align Layers");
    gd.addChoice("mode :", modeStrings, modeStrings[LINEAR]);
    gd.addMessage("Layer Range:");
    final int sel = l.getParent().indexOf(l);
    gd.addChoice("first :", layerTitles, layerTitles[sel]);
    gd.addChoice("reference :", noneAndLayerTitles, layerTitles[sel]);
    gd.addChoice("last :", layerTitles, layerTitles[sel]);
    gd.addStringField("Use only images whose title matches:", "", 30);
    gd.addCheckbox("Use visible images only", true);
    gd.addCheckbox("propagate transform to first layer", propagateTransformBefore);
    gd.addCheckbox("propagate transform to last layer", propagateTransformAfter);
    final Vector<?> v = gd.getChoices();
    final Choice cstart = (Choice) v.get(v.size() - 3);
    final Choice cref = (Choice) v.get(v.size() - 2);
    final Choice cend = (Choice) v.get(v.size() - 1);
    final ItemListener startEndListener = new ItemListener() {

        @Override
        public void itemStateChanged(final ItemEvent ie) {
            final int startIndex = cstart.getSelectedIndex();
            final int endIndex = cend.getSelectedIndex();
            final int refIndex = cref.getSelectedIndex();
            final int min = Math.min(startIndex, endIndex);
            final int max = Math.max(startIndex, endIndex);
            if (cref.getSelectedIndex() > 0) {
                cref.select(Math.min(max + 1, Math.max(min + 1, refIndex)));
            }
        }
    };
    cstart.addItemListener(startEndListener);
    cend.addItemListener(startEndListener);
    cref.addItemListener(new ItemListener() {

        @Override
        public void itemStateChanged(final ItemEvent ie) {
            final int startIndex = cstart.getSelectedIndex();
            final int endIndex = cend.getSelectedIndex();
            final int refIndex = cref.getSelectedIndex() - 1;
            final int min = Math.min(startIndex, endIndex);
            final int max = Math.max(startIndex, endIndex);
            if (refIndex >= 0) {
                if (refIndex < min) {
                    if (startIndex <= endIndex)
                        cstart.select(refIndex);
                    else
                        cend.select(refIndex);
                } else if (refIndex > max) {
                    if (startIndex <= endIndex)
                        cend.select(refIndex);
                    else
                        cstart.select(refIndex);
                }
            }
        }
    });
    gd.showDialog();
    if (gd.wasCanceled())
        return;
    mode = gd.getNextChoiceIndex();
    final int first = gd.getNextChoiceIndex();
    final int ref = gd.getNextChoiceIndex() - 1;
    final int last = gd.getNextChoiceIndex();
    final String toMatch1 = gd.getNextString().trim();
    final String toMatch2 = 0 == toMatch1.length() ? null : ".*" + toMatch1 + ".*";
    final boolean visibleOnly = gd.getNextBoolean();
    propagateTransformBefore = gd.getNextBoolean();
    propagateTransformAfter = gd.getNextBoolean();
    final Filter<Patch> filter = new Filter<Patch>() {

        @Override
        public final boolean accept(final Patch patch) {
            if (visibleOnly && !patch.isVisible())
                return false;
            if (null != toMatch2 && !patch.getTitle().matches(toMatch2))
                return false;
            return true;
        }
    };
    if (mode == ELASTIC)
        new ElasticLayerAlignment().exec(l.getParent(), first, last, ref, propagateTransformBefore, propagateTransformAfter, fov, filter);
    else if (mode == LINEAR)
        new RegularizedAffineLayerAlignment().exec(l.getParent(), first, last, ref, propagateTransformBefore, propagateTransformAfter, fov, filter);
    else {
        final GenericDialog gd2 = new GenericDialog("Align Layers");
        Align.param.addFields(gd2);
        gd2.showDialog();
        if (gd2.wasCanceled())
            return;
        Align.param.readFields(gd2);
        if (mode == BUNWARPJ && !elasticParam.showDialog())
            return;
        // From ref to first:
        if (ref - first > 0) {
            if (mode == BUNWARPJ)
                alignLayersNonLinearlyJob(l.getParent(), ref, first, propagateTransformBefore, fov, filter);
            else
                alignLayersLinearlyJob(l.getParent(), ref, first, propagateTransformBefore, fov, filter);
        }
        // From ref to last:
        if (last - ref > 0) {
            if (mode == BUNWARPJ)
                alignLayersNonLinearlyJob(l.getParent(), ref, last, propagateTransformAfter, fov, filter);
            else
                alignLayersLinearlyJob(l.getParent(), ref, last, propagateTransformAfter, fov, filter);
        }
    }
}

Example 10 with Filter

use of ini.trakem2.utils.Filter in project TrakEM2 by trakem2.

the class MatchIntensities method run.

/**
 * @param layers
 * @param radius
 * @param scale
 * @param numCoefficients
 * @param lambda1
 * @param lambda2
 * @param neighborWeight
 * @param roi
 */
public <M extends Model<M> & Affine1D<M>> void run(final List<Layer> layers, final int radius, final double scale, final int numCoefficients, final double lambda1, final double lambda2, final double neighborWeight, final Rectangle roi) throws InterruptedException, ExecutionException {
    final int firstLayerIndex = layerset.getLayerIndex(layers.get(0).getId());
    final int lastLayerIndex = layerset.getLayerIndex(layers.get(layers.size() - 1).getId());
    // final PointMatchFilter filter = new RansacRegressionFilter();
    final PointMatchFilter filter = new RansacRegressionReduceFilter();
    /* collect patches */
    Utils.log("Collecting patches ... ");
    final ArrayList<Patch> patches = new ArrayList<Patch>();
    for (final Layer layer : layers) patches.addAll((Collection) layer.getDisplayables(Patch.class, roi));
    /* delete existing intensity coefficients */
    Utils.log("Clearing existing intensity maps ... ");
    for (final Patch p : patches) p.clearIntensityMap();
    /* generate coefficient tiles for all patches
		 * TODO consider offering alternative models */
    final HashMap<Patch, ArrayList<Tile<? extends M>>> coefficientsTiles = (HashMap) generateCoefficientsTiles(patches, new InterpolatedAffineModel1D<InterpolatedAffineModel1D<AffineModel1D, TranslationModel1D>, IdentityModel>(new InterpolatedAffineModel1D<AffineModel1D, TranslationModel1D>(new AffineModel1D(), new TranslationModel1D(), lambda1), new IdentityModel(), lambda2), numCoefficients * numCoefficients);
    /* completed patches */
    final HashSet<Patch> completedPatches = new HashSet<Patch>();
    /* collect patch pairs */
    Utils.log("Collecting patch pairs ... ");
    final ArrayList<ValuePair<Patch, Patch>> patchPairs = new ArrayList<ValuePair<Patch, Patch>>();
    for (final Patch p1 : patches) {
        completedPatches.add(p1);
        final Rectangle box1 = p1.getBoundingBox().intersection(roi);
        final ArrayList<Patch> p2s = new ArrayList<Patch>();
        /* across adjacent layers */
        final int layerIndex = layerset.getLayerIndex(p1.getLayer().getId());
        for (int i = Math.max(firstLayerIndex, layerIndex - radius); i <= Math.min(lastLayerIndex, layerIndex + radius); ++i) {
            final Layer layer = layerset.getLayer(i);
            if (layer != null)
                p2s.addAll((Collection) layer.getDisplayables(Patch.class, box1));
        }
        for (final Patch p2 : p2s) {
            /*
				 * if this patch had been processed earlier, all matches are
				 * already in
				 */
            if (completedPatches.contains(p2))
                continue;
            patchPairs.add(new ValuePair<Patch, Patch>(p1, p2));
        }
    }
    final int numThreads = Integer.parseInt(layerset.getProperty("n_mipmap_threads", Integer.toString(Runtime.getRuntime().availableProcessors())));
    Utils.log("Matching intensities using " + numThreads + " threads ... ");
    final ExecutorService exec = Executors.newFixedThreadPool(numThreads);
    final ArrayList<Future<?>> futures = new ArrayList<Future<?>>();
    for (final ValuePair<Patch, Patch> patchPair : patchPairs) {
        futures.add(exec.submit(new Matcher(roi, patchPair, (HashMap) coefficientsTiles, filter, scale, numCoefficients)));
    }
    for (final Future<?> future : futures) future.get();
    /* connect tiles within patches */
    Utils.log("Connecting coefficient tiles in the same patch  ... ");
    for (final Patch p1 : completedPatches) {
        /* get the coefficient tiles */
        final ArrayList<Tile<? extends M>> p1CoefficientsTiles = coefficientsTiles.get(p1);
        for (int y = 1; y < numCoefficients; ++y) {
            final int yr = numCoefficients * y;
            final int yr1 = yr - numCoefficients;
            for (int x = 0; x < numCoefficients; ++x) {
                identityConnect(p1CoefficientsTiles.get(yr1 + x), p1CoefficientsTiles.get(yr + x), neighborWeight);
            }
        }
        for (int y = 0; y < numCoefficients; ++y) {
            final int yr = numCoefficients * y;
            for (int x = 1; x < numCoefficients; ++x) {
                final int yrx = yr + x;
                identityConnect(p1CoefficientsTiles.get(yrx), p1CoefficientsTiles.get(yrx - 1), neighborWeight);
            }
        }
    }
    /* optimize */
    Utils.log("Optimizing ... ");
    final TileConfiguration tc = new TileConfiguration();
    for (final ArrayList<Tile<? extends M>> coefficients : coefficientsTiles.values()) {
        // for ( final Tile< ? > t : coefficients )
        // if ( t.getMatches().size() == 0 )
        // IJ.log( "bang" );
        tc.addTiles(coefficients);
    }
    try {
        tc.optimize(0.01f, iterations, iterations, 0.75f);
    } catch (final NotEnoughDataPointsException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    } catch (final IllDefinedDataPointsException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
    /* save coefficients */
    final double[] ab = new double[2];
    final FSLoader loader = (FSLoader) layerset.getProject().getLoader();
    final String itsDir = loader.getUNUIdFolder() + "trakem2.its/";
    for (final Entry<Patch, ArrayList<Tile<? extends M>>> entry : coefficientsTiles.entrySet()) {
        final FloatProcessor as = new FloatProcessor(numCoefficients, numCoefficients);
        final FloatProcessor bs = new FloatProcessor(numCoefficients, numCoefficients);
        final Patch p = entry.getKey();
        final double min = p.getMin();
        final double max = p.getMax();
        final ArrayList<Tile<? extends M>> tiles = entry.getValue();
        for (int i = 0; i < numCoefficients * numCoefficients; ++i) {
            final Tile<? extends M> t = tiles.get(i);
            final Affine1D<?> affine = t.getModel();
            affine.toArray(ab);
            /* coefficients mapping into existing [min, max] */
            as.setf(i, (float) ab[0]);
            bs.setf(i, (float) ((max - min) * ab[1] + min - ab[0] * min));
        }
        final ImageStack coefficientsStack = new ImageStack(numCoefficients, numCoefficients);
        coefficientsStack.addSlice(as);
        coefficientsStack.addSlice(bs);
        final String itsPath = itsDir + FSLoader.createIdPath(Long.toString(p.getId()), "it", ".tif");
        new File(itsPath).getParentFile().mkdirs();
        IJ.saveAs(new ImagePlus("", coefficientsStack), "tif", itsPath);
    }
    /* update mipmaps */
    for (final Patch p : patches) p.getProject().getLoader().decacheImagePlus(p.getId());
    final ArrayList<Future<Boolean>> mipmapFutures = new ArrayList<Future<Boolean>>();
    for (final Patch p : patches) mipmapFutures.add(p.updateMipMaps());
    for (final Future<Boolean> f : mipmapFutures) f.get();
    Utils.log("Matching intensities done.");
}