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

Example 1 with CircularSequence

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