Examples with VectorString2D ini.trakem2.vector.VectorString2D used on opensource projects

Example 1 with VectorString2D

use of ini.trakem2.vector.VectorString2D in project TrakEM2 by trakem2.

the class Profile method getPerimeter2D.

private VectorString2D getPerimeter2D(final Calibration cal) {
    if (-1 == n_points)
        setupForDisplay();
    if (0 == n_points)
        return null;
    if (0 == p_i[0].length)
        generateInterpolatedPoints(0.05);
    final double[][] pi = transformPoints(p_i);
    VectorString2D sv = null;
    try {
        sv = new VectorString2D(pi[0], pi[1], this.layer.getZ(), this.closed);
    } catch (final Exception e) {
        IJError.print(e);
    }
    if (null != cal)
        sv.calibrate(cal);
    return sv;
}

Example 2 with VectorString2D

use of ini.trakem2.vector.VectorString2D in project TrakEM2 by trakem2.

the class AreaUtils method singularInterpolation.

/**
 * Extract the Area of the image for the given pixel value.
 *  ThresholdToSelection is way faster than this, just use it.
 *  It's BROKEN do not use.
 */
/*
	static public final Area extractArea(final ImageProcessor ip, final float val) {
		final int height = ip.getHeight();
		final int width = ip.getWidth();

		final Area area = new Area();
		final ArrayList<Area> segments = new ArrayList<Area>();
		final Rectangle box = new Rectangle(0, 0, 1, 1);

		for (int y=0; y<height; y++) {
			float prev = ip.getPixelValue(0, y);
			box.x = 0;
			box.y = y;
			box.width = val == prev ? 1 : 0;

			for (int x=1; x<width; x++) {

				float pix = ip.getPixelValue(x, y);

				if (val == pix) {
					if (pix == prev) {
						box.width++;
						continue;
					}
					// Else, add previous one
					segments.add(new Area(box));
					// ... and start a new box
					box.x = x;
					box.y = y;
					box.width = 1;
					prev = pix;
				}
			}

			// At end of line, add the last
			segments.add(new Area(box));

			// Join a few
			if (segments.size() > 32) {
				final Area a = new Area(segments.get(0));
				final int len = segments.size();
				for (int i=1; i<len; i++) a.add(segments.get(i));
				area.add(a);
				segments.clear();
			}
		}

		if (segments.size() > 0) {
			final Area a = new Area(segments.get(0));
			final int len = segments.size();
			for (int i=1; i<len; i++) a.add(segments.get(i));
			area.add(a);
		}

		return area;
	}
	*/
/**
 * Interpolate areas only if they are made of a single shape each.
 *  Assumes that areas are in the same coordinate system.
 * @throws Exception
 */
public static final Area[] singularInterpolation(final Area a1, final Area a2, final int nInterpolates) throws Exception {
    if (!a1.isSingular() || !a2.isSingular()) {
        return null;
    }
    final VectorString2D vs1 = M.asVectorString2D(M.getPolygons(a1).iterator().next(), 0);
    final VectorString2D vs2 = M.asVectorString2D(M.getPolygons(a2).iterator().next(), 1);
    final Editions ed = new Editions(vs1, vs2, Math.min(vs1.getAverageDelta(), vs2.getAverageDelta()), true);
    final double[][][] d = SkinMaker.getMorphedPerimeters(vs1, vs2, nInterpolates, ed);
    final Area[] a = new Area[d.length];
    for (int i = 0; i < d.length; i++) {
        final double[] x = d[i][0];
        final double[] y = d[i][1];
        final int[] xi = new int[x.length];
        final int[] yi = new int[y.length];
        for (int k = 0; k < x.length; k++) {
            xi[k] = (int) x[k];
            yi[k] = (int) y[k];
        }
        a[i] = new Area(new Polygon(xi, yi, xi.length));
    }
    return a;
}

Example 3 with VectorString2D

use of ini.trakem2.vector.VectorString2D in project TrakEM2 by trakem2.

the class AreaList method measure.

/**
 * Returns a double array with 0=volume, 1=lower_bound_surface, 2=upper_bound_surface_smoothed, 3=upper_bound_surface, 4=max_diameter, 5=all_tops_and_bottoms
 *  All measures are approximate.
 *  [0] Volume: sum(area * thickness) for all sections
 *  [1] Lower Bound Surface: measure area per section, compute radius of circumference of identical area, compute then area of the sides of the truncated cone of height thickness, for each section. Plus top and bottom areas when visiting sections without a painted area.
 *  [2] Upper Bound Surface Smoothed: measure smoothed perimeter lengths per section, multiply by thickness to get lateral area. Plus tops and bottoms.
 *  [3] Upper Bound Surface: measure raw pixelated perimeter lengths per section, multiply by thickness to get lateral area. Plus top and bottoms.
 *  [4] Maximum diameter: longest distance between any two points in the contours of all painted areas.
 *  [5] All tops and bottoms: Sum of all included surface areas that are not part of side area.
 *  [6] X coordinate of the center of mass.
 *  [7] Y coordinate of the center of mass.
 *  [8] Z coordinate of the center of mass.
 */
public double[] measure() {
    // zeros
    if (0 == ht_areas.size())
        return new double[6];
    // prepare suitable transform
    final AffineTransform aff = new AffineTransform(this.at);
    AffineTransform aff2 = new AffineTransform();
    // remove translation (for no reason other than historical, and that it may
    // help avoid numerical overflows)
    final Rectangle box = getBoundingBox(null);
    aff2.translate(-box.x, -box.y);
    aff.preConcatenate(aff2);
    aff2 = null;
    double volume = 0;
    double lower_bound_surface_h = 0;
    double upper_bound_surface = 0;
    double upper_bound_surface_smoothed = 0;
    double prev_surface = 0;
    double prev_perimeter = 0;
    double prev_smooth_perimeter = 0;
    double prev_thickness = 0;
    // i.e. surface area that is not part of the side area
    double all_tops_and_bottoms = 0;
    final Calibration cal = layer_set.getCalibration();
    final double pixelWidth = cal.pixelWidth;
    final double pixelHeight = cal.pixelHeight;
    // Put areas in order of their layer index:
    final TreeMap<Integer, Area> ias = new TreeMap<Integer, Area>();
    for (final Map.Entry<Long, Area> e : ht_areas.entrySet()) {
        final int ilayer = layer_set.indexOf(layer_set.getLayer(e.getKey()));
        if (-1 == ilayer) {
            Utils.log("Could not find a layer with id " + e.getKey());
            continue;
        }
        ias.put(ilayer, e.getValue());
    }
    final ArrayList<Layer> layers = layer_set.getLayers();
    int last_layer_index = -1;
    final ArrayList<Point3f> points = new ArrayList<Point3f>();
    final float[] coords = new float[6];
    final float fpixelWidth = (float) pixelWidth;
    final float fpixelHeight = (float) pixelHeight;
    final float resampling_delta = project.getProperty("measurement_resampling_delta", 1.0f);
    // for each area, measure its area and its perimeter, to compute volume and surface
    for (final Map.Entry<Integer, Area> e : ias.entrySet()) {
        // fetch Layer
        final int layer_index = e.getKey();
        if (layer_index > layers.size()) {
            Utils.log("Could not find a layer at index " + layer_index);
            continue;
        }
        final Layer la = layers.get(layer_index);
        // fetch Area
        Area area = e.getValue();
        if (UNLOADED == area)
            area = loadLayer(la.getId());
        // Transform area to world coordinates
        area = area.createTransformedArea(aff);
        // measure surface
        final double pixel_area = Math.abs(AreaCalculations.area(area.getPathIterator(null)));
        final double surface = pixel_area * pixelWidth * pixelHeight;
        // Utils.log2(layer_index + " pixel_area: " + pixel_area + "  surface " + surface);
        // measure volume
        // the last one is NOT pixelDepth because layer thickness and Z are in pixels
        final double thickness = la.getThickness() * pixelWidth;
        volume += surface * thickness;
        final double pix_perimeter = AreaCalculations.circumference(area.getPathIterator(null));
        final double perimeter = pix_perimeter * pixelWidth;
        double smooth_perimeter = 0;
        // smoothed perimeter:
        {
            double smooth_pix_perimeter = 0;
            for (final Polygon pol : M.getPolygons(area)) {
                try {
                    if (pol.npoints < 5) {
                        // No point in smoothing out such a short polygon:
                        // (Plus can't convolve it with a gaussian that needs 5 points adjacent)
                        smooth_perimeter += new PolygonRoi(pol, PolygonRoi.POLYGON).getLength();
                        continue;
                    }
                    /*
						// Works but it is not the best smoothing of the Area's countour
						double[] xp = new double[pol.npoints];
						double[] yp = new double[pol.npoints];
						for (int p=0; p<pol.npoints; p++) {
							xp[p] = pol.xpoints[p];
							yp[p] = pol.ypoints[p];
						}
						VectorString2D v = new VectorString2D(xp, yp, 0, true);
						v.resample(resampling_delta);
						smooth_pix_perimeter += v.length() * resampling_delta;
						*/
                    // The best solution I've found:
                    // 1. Run getInterpolatedPolygon with an interval of 1 to get a point at every pixel
                    // 2. convolve with a gaussian
                    // Resample to 1 so that at every one pixel of the contour there is a point
                    FloatPolygon fpol = new FloatPolygon(new float[pol.npoints], new float[pol.npoints], pol.npoints);
                    for (int i = 0; i < pol.npoints; ++i) {
                        fpol.xpoints[i] = pol.xpoints[i];
                        fpol.ypoints[i] = pol.ypoints[i];
                    }
                    fpol = M.createInterpolatedPolygon(fpol, 1, false);
                    final FloatPolygon fp;
                    if (fpol.npoints < 5) {
                        smooth_pix_perimeter += fpol.getLength(false);
                        fp = fpol;
                    } else {
                        // Convolve with a sigma of 1 to smooth it out
                        final FloatPolygon gpol = new FloatPolygon(new float[fpol.npoints], new float[fpol.npoints], fpol.npoints);
                        final CircularSequence seq = new CircularSequence(fpol.npoints);
                        M.convolveGaussianSigma1(fpol.xpoints, gpol.xpoints, seq);
                        M.convolveGaussianSigma1(fpol.ypoints, gpol.ypoints, seq);
                        // Resample it to the desired resolution (also facilitates measurement: npoints * resampling_delta)
                        if (gpol.npoints > resampling_delta) {
                            fp = M.createInterpolatedPolygon(gpol, resampling_delta, false);
                        } else {
                            fp = gpol;
                        }
                        // Measure perimeter: last line segment is potentially shorter or longer than resampling_delta
                        smooth_pix_perimeter += (fp.npoints - 1) * resampling_delta + Math.sqrt(Math.pow(fp.xpoints[0] - fp.xpoints[fp.npoints - 1], 2) + Math.pow(fp.ypoints[0] - fp.ypoints[fp.npoints - 1], 2));
                    }
                // TEST:
                // ij.plugin.frame.RoiManager.getInstance().addRoi(new PolygonRoi(fp, PolygonRoi.POLYGON));
                // TESTING: make a polygon roi and show it
                // ... just in case to see that resampling works as expected, without weird endings
                /*
						int[] x = new int[v.length()];
						int[] y = new int[x.length];
						double[] xd = v.getPoints(0);
						double[] yd = v.getPoints(1);
						for (int p=0; p<x.length; p++) {
							x[p] = (int)xd[p];
							y[p] = (int)yd[p];
						}
						PolygonRoi proi = new PolygonRoi(x, y, x.length, PolygonRoi.POLYGON);
						Rectangle b = proi.getBounds();
						for (int p=0; p<x.length; p++) {
							x[p] -= b.x;
							y[p] -= b.y;
						}
						ImagePlus imp = new ImagePlus("test", new ByteProcessor(b.width, b.height));
						imp.setRoi(new PolygonRoi(x, y, x.length, PolygonRoi.POLYGON));
						imp.show();
						*/
                } catch (final Exception le) {
                    le.printStackTrace();
                }
            }
            smooth_perimeter = smooth_pix_perimeter * pixelWidth;
        }
        if (-1 == last_layer_index) {
            // Start of the very first continuous set:
            lower_bound_surface_h += surface;
            upper_bound_surface += surface;
            upper_bound_surface_smoothed += surface;
            all_tops_and_bottoms += surface;
        } else if (layer_index - last_layer_index > 1) {
            // End of a continuous set ...
            // Sum the last surface and its side:
            // (2x + 2x) / 2   ==   2x
            lower_bound_surface_h += prev_surface + prev_thickness * 2 * Math.sqrt(prev_surface * Math.PI);
            upper_bound_surface += prev_surface + prev_perimeter * prev_thickness;
            upper_bound_surface_smoothed += prev_surface + prev_smooth_perimeter * prev_thickness;
            all_tops_and_bottoms += prev_surface;
            // ... and start of a new set
            lower_bound_surface_h += surface;
            upper_bound_surface += surface;
            upper_bound_surface_smoothed += surface;
            all_tops_and_bottoms += surface;
        } else {
            // Continuation of a set: use this Area and the previous as continuous
            final double diff_surface = Math.abs(prev_surface - surface);
            upper_bound_surface += prev_perimeter * (prev_thickness / 2) + perimeter * (prev_thickness / 2) + diff_surface;
            upper_bound_surface_smoothed += prev_smooth_perimeter * (prev_thickness / 2) + smooth_perimeter * (prev_thickness / 2) + diff_surface;
            // Compute area of the mantle of the truncated cone defined by the radiuses of the circles of same area as the two areas
            // PI * s * (r1 + r2) where s is the hypothenusa
            final double r1 = Math.sqrt(prev_surface / Math.PI);
            final double r2 = Math.sqrt(surface / Math.PI);
            final double hypothenusa = Math.sqrt(Math.pow(Math.abs(r1 - r2), 2) + Math.pow(thickness, 2));
            lower_bound_surface_h += Math.PI * hypothenusa * (r1 + r2);
            // Adjust volume too:
            volume += diff_surface * prev_thickness / 2;
        }
        // store for next iteration:
        prev_surface = surface;
        prev_perimeter = perimeter;
        prev_smooth_perimeter = smooth_perimeter;
        last_layer_index = layer_index;
        prev_thickness = thickness;
        // Iterate points:
        final float z = (float) la.getZ();
        for (final PathIterator pit = area.getPathIterator(null); !pit.isDone(); pit.next()) {
            switch(pit.currentSegment(coords)) {
                case PathIterator.SEG_MOVETO:
                case PathIterator.SEG_LINETO:
                case PathIterator.SEG_CLOSE:
                    points.add(new Point3f(coords[0] * fpixelWidth, coords[1] * fpixelHeight, z * fpixelWidth));
                    break;
                default:
                    Utils.log2("WARNING: unhandled seg type.");
                    break;
            }
        }
    }
    // finish last:
    lower_bound_surface_h += prev_surface + prev_perimeter * prev_thickness;
    upper_bound_surface += prev_surface + prev_perimeter * prev_thickness;
    upper_bound_surface_smoothed += prev_surface + prev_smooth_perimeter * prev_thickness;
    all_tops_and_bottoms += prev_surface;
    // Compute maximum diameter
    final boolean measure_largest_diameter = project.getBooleanProperty("measure_largest_diameter");
    double max_diameter_sq = measure_largest_diameter ? 0 : Double.NaN;
    final int lp = points.size();
    final Point3f c;
    if (lp > 0) {
        // center of mass
        c = new Point3f(points.get(0));
        for (int i = 0; i < lp; i++) {
            final Point3f p = points.get(i);
            if (measure_largest_diameter) {
                for (int j = i; j < lp; j++) {
                    final double len = p.distanceSquared(points.get(j));
                    if (len > max_diameter_sq)
                        max_diameter_sq = len;
                }
            }
            if (0 == i)
                continue;
            c.x += p.x;
            c.y += p.y;
            c.z += p.z;
        }
    } else {
        c = new Point3f(Float.NaN, Float.NaN, Float.NaN);
    }
    // Translate the center of mass
    c.x = box.x + c.x / lp;
    c.y = box.y + c.y / lp;
    c.z /= lp;
    return new double[] { volume, lower_bound_surface_h, upper_bound_surface_smoothed, upper_bound_surface, Math.sqrt(max_diameter_sq), all_tops_and_bottoms, c.x, c.y, c.z };
}

Example 4 with VectorString2D

use of ini.trakem2.vector.VectorString2D in project TrakEM2 by trakem2.

the class Profile method makeTriangles.

/**
 * Make a mesh as a calibrated list of 3D triangles.
 */
private static List<Point3f> makeTriangles(final Profile[] p, final double scale) {
    try {
        final VectorString2D[] sv = new VectorString2D[p.length];
        // dummy initialization
        boolean closed = true;
        final Calibration cal = p[0].getLayerSet().getCalibrationCopy();
        cal.pixelWidth *= scale;
        cal.pixelHeight *= scale;
        for (int i = 0; i < p.length; i++) {
            if (0 == p[i].n_points)
                continue;
            if (0 == i)
                closed = p[i].closed;
            else if (p[i].closed != closed) {
                Utils.log2("All profiles should be either open or closed, not mixed.");
                return null;
            }
            sv[i] = p[i].getPerimeter2D(cal);
        }
        return SkinMaker.generateTriangles(sv, -1, -1, closed);
    } catch (final Exception e) {
        IJError.print(e);
    }
    return null;
}