use of mpicbg.models.ErrorStatistic in project TrakEM2 by trakem2.
the class StitchingTEM method optimizeTileConfiguration.
/**
* Optimize tile configuration by removing bad matches
*
* @param tiles complete list of tiles
* @param fixed_tiles list of fixed tiles
* @param param phase correlation parameters
*/
public static void optimizeTileConfiguration(final ArrayList<AbstractAffineTile2D<?>> tiles, final ArrayList<AbstractAffineTile2D<?>> fixed_tiles, final PhaseCorrelationParam param) {
// Run optimization
if (fixed_tiles.isEmpty())
fixed_tiles.add(tiles.get(0));
// with default parameters
boolean proceed = true;
while (proceed) {
Align.optimizeTileConfiguration(new Align.ParamOptimize(), tiles, fixed_tiles);
/* get all transfer errors */
final ErrorStatistic e = new ErrorStatistic(tiles.size() + 1);
for (final AbstractAffineTile2D<?> t : tiles) t.update();
for (final AbstractAffineTile2D<?> t : tiles) {
for (final PointMatch p : t.getMatches()) {
e.add(p.getDistance());
}
}
/* remove the worst if there is one */
if (e.max > param.mean_factor * e.mean) {
A: for (final AbstractAffineTile2D<?> t : tiles) {
for (final PointMatch p : t.getMatches()) {
if (p.getDistance() >= e.max) {
final Tile<?> o = t.findConnectedTile(p);
t.removeConnectedTile(o);
o.removeConnectedTile(t);
// Utils.log2( "Removing bad match from configuration, error = " + e.max );
break A;
}
}
}
} else
proceed = false;
}
}
use of mpicbg.models.ErrorStatistic 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.");
}
use of mpicbg.models.ErrorStatistic in project TrakEM2 by trakem2.
the class BlockMatchPairCallable method call.
@Override
public BlockMatchResults call() throws Exception {
final ArrayList<PointMatch> pm12 = new ArrayList<PointMatch>();
final ArrayList<PointMatch> pm21 = new ArrayList<PointMatch>();
System.out.println("BMC rev 0: " + pair.a + " " + pair.b);
final Pair<FloatProcessor, FloatProcessor> pair1 = makeFlatImage(layer1, AlignmentUtils.filterPatches(layer1, filter), box, param.layerScale);
final Pair<FloatProcessor, FloatProcessor> pair2 = makeFlatImage(layer2, AlignmentUtils.filterPatches(layer2, filter), box, param.layerScale);
final FloatProcessor ip1 = pair1.a;
final FloatProcessor ip1Mask = pair1.b;
final FloatProcessor ip2 = pair2.a;
final FloatProcessor ip2Mask = pair2.b;
final AbstractModel<?> localSmoothnessFilterModel = Util.createModel(param.localModelIndex);
final int blockRadius = Math.max(16, mpicbg.util.Util.roundPos(param.layerScale * param.blockRadius));
/* scale pixel distances */
final int searchRadius = (int) Math.round(param.layerScale * param.searchRadius);
final double localRegionSigma = param.layerScale * param.localRegionSigma;
final double maxLocalEpsilon = param.layerScale * param.maxLocalEpsilon;
if (!layer1Fixed) {
BlockMatching.matchByMaximalPMCC(ip1, ip2, ip1Mask, ip2Mask, 1.0, ((InvertibleCoordinateTransform) pair.c).createInverse(), blockRadius, blockRadius, searchRadius, searchRadius, param.minR, param.rodR, param.maxCurvatureR, v1, pm12, new ErrorStatistic(1));
if (Thread.interrupted()) {
throw new InterruptedException("Block matching interrupted.");
}
if (param.useLocalSmoothnessFilter) {
localSmoothnessFilterModel.localSmoothnessFilter(pm12, pm12, localRegionSigma, maxLocalEpsilon, param.maxLocalTrust);
}
}
if (!layer2Fixed) {
BlockMatching.matchByMaximalPMCC(ip2, ip1, ip2Mask, ip1Mask, 1.0f, pair.c, blockRadius, blockRadius, searchRadius, searchRadius, param.minR, param.rodR, param.maxCurvatureR, v2, pm21, new ErrorStatistic(1));
if (Thread.interrupted()) {
throw new InterruptedException("Block matching interrupted.");
}
if (param.useLocalSmoothnessFilter) {
localSmoothnessFilterModel.localSmoothnessFilter(pm21, pm21, localRegionSigma, maxLocalEpsilon, param.maxLocalTrust);
}
}
return new BlockMatchResults(v1, v2, pm12, pm21, layer1Fixed, layer2Fixed, pair);
}
Aggregations