Search in sources :

Example 16 with Coordinate

use of ini.trakem2.display.Coordinate in project TrakEM2 by trakem2.

the class Project method saveWithoutCoordinateTransforms.

/**
 * Save an XML file that is stripped of coordinate transforms,
 * and merely refers to them by the 'ct_id' attribute of each 't2_patch' element;
 * this method will NOT overwrite the XML file but save into a new one,
 * which is chosen from a file dialog.
 */
public String saveWithoutCoordinateTransforms() {
    XMLOptions options = new XMLOptions();
    options.overwriteXMLFile = false;
    options.export_images = false;
    options.include_coordinate_transform = false;
    options.patches_dir = null;
    return loader.saveAs(this, options);
}
Also used : XMLOptions(ini.trakem2.persistence.XMLOptions)

Example 17 with Coordinate

use of ini.trakem2.display.Coordinate 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 };
}
Also used : PathIterator(java.awt.geom.PathIterator) Rectangle(java.awt.Rectangle) ArrayList(java.util.ArrayList) PolygonRoi(ij.gui.PolygonRoi) Point3f(org.scijava.vecmath.Point3f) FloatPolygon(ij.process.FloatPolygon) Polygon(java.awt.Polygon) Calibration(ij.measure.Calibration) TreeMap(java.util.TreeMap) Point(java.awt.Point) USHORTPaint(ini.trakem2.display.paint.USHORTPaint) NoninvertibleTransformException(java.awt.geom.NoninvertibleTransformException) Area(java.awt.geom.Area) AffineTransform(java.awt.geom.AffineTransform) Map(java.util.Map) HashMap(java.util.HashMap) TreeMap(java.util.TreeMap) FloatPolygon(ij.process.FloatPolygon) CircularSequence(ini.trakem2.utils.CircularSequence)

Example 18 with Coordinate

use of ini.trakem2.display.Coordinate in project TrakEM2 by trakem2.

the class ProjectThing method fixZOrdering.

/**
 * If this is of type profile_list, order all children Profile by their layer's Z coordinate, ascending.
 */
public boolean fixZOrdering() {
    if (!this.template.getType().equals("profile_list"))
        return false;
    // no need
    if (null == al_children || al_children.size() < 2)
        return true;
    // fix Z ordering in the tree
    final HashMap<Double, ArrayList<ProjectThing>> ht = new HashMap<Double, ArrayList<ProjectThing>>();
    synchronized (al_children) {
        for (ProjectThing child : al_children) {
            Profile p = (Profile) child.object;
            Layer layer = p.getLayer();
            // contortions: there could be more than one profile in the same layer
            Double z = new Double(layer.getZ());
            ArrayList<ProjectThing> al = ht.get(z);
            if (null == al) {
                al = new ArrayList<ProjectThing>();
                al.add(child);
                ht.put(z, al);
            } else {
                al.add(child);
            }
        }
        Double[] zs = new Double[ht.size()];
        ht.keySet().toArray(zs);
        Arrays.sort(zs);
        al_children.clear();
        for (int i = 0; i < zs.length; i++) {
            al_children.addAll(ht.get(zs[i]));
        }
    }
    return true;
}
Also used : HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) Layer(ini.trakem2.display.Layer) Profile(ini.trakem2.display.Profile)

Example 19 with Coordinate

use of ini.trakem2.display.Coordinate in project TrakEM2 by trakem2.

the class Search method executeSearch.

private void executeSearch() {
    final Project project = Project.getProjects().get(projects.getSelectedIndex());
    if (null == project) {
        // Should not happen
        return;
    }
    Bureaucrat.createAndStart(new Worker.Task("Searching") {

        public void exec() {
            String pattern = search_field.getText();
            if (null == pattern || 0 == pattern.length()) {
                return;
            }
            // fix pattern
            final String typed_pattern = pattern;
            // I hate java
            final StringBuilder sb = new StringBuilder();
            if (!pattern.startsWith("^"))
                sb.append("^.*");
            for (int i = 0; i < pattern.length(); i++) {
                sb.append(pattern.charAt(i));
            }
            if (!pattern.endsWith("$"))
                sb.append(".*$");
            pattern = sb.toString();
            final Pattern pat = Pattern.compile(pattern, Pattern.CASE_INSENSITIVE | Pattern.MULTILINE | Pattern.DOTALL);
            // Utils.log2("pattern after: " + pattern);
            final ArrayList<DBObject> al = new ArrayList<DBObject>();
            // Utils.log("types[pulldown] = " +
            // types[pulldown.getSelectedIndex()]);
            find(project.getRootLayerSet(), al, types[pulldown.getSelectedIndex()]);
            // Utils.log2("found labels: " + al.size());
            if (0 == al.size())
                return;
            final Vector<DBObject> v_obs = new Vector<DBObject>();
            final Vector<String> v_txt = new Vector<String>();
            final Vector<Coordinate<?>> v_co = new Vector<Coordinate<?>>();
            Coordinate<?> co = null;
            for (final DBObject dbo : al) {
                if (Thread.currentThread().isInterrupted()) {
                    return;
                }
                boolean matched = false;
                // Search in its title
                Displayable d = null;
                if (dbo instanceof Displayable) {
                    d = (Displayable) dbo;
                }
                String txt;
                String meaningful_title = null;
                if (null == d || Patch.class == d.getClass())
                    txt = dbo.getTitle();
                else {
                    txt = meaningful_title = dbo.getProject().getMeaningfulTitle(d);
                }
                if (null == txt || 0 == txt.trim().length())
                    continue;
                matched = pat.matcher(txt).matches();
                if (!matched && null != d) {
                    // Search also in its annotation
                    txt = d.getAnnotation();
                    if (null != txt)
                        matched = pat.matcher(txt).matches();
                }
                if (!matched) {
                    // Search also in its toString()
                    txt = dbo.toString();
                    matched = pat.matcher(txt).matches();
                }
                if (!matched) {
                    // Search also in its id
                    txt = Long.toString(dbo.getId());
                    matched = pat.matcher(txt).matches();
                    if (matched)
                        txt = "id: #" + txt;
                }
                if (!matched && null != d) {
                    // Search also in its properties
                    Map<String, String> props = d.getProperties();
                    if (null != props) {
                        for (final Map.Entry<String, String> e : props.entrySet()) {
                            if (pat.matcher(e.getKey()).matches() || pat.matcher(e.getValue()).matches()) {
                                matched = true;
                                txt = e.getKey() + " => " + e.getValue() + " [property]";
                                break;
                            }
                        }
                    }
                    if (!matched) {
                        Map<Displayable, Map<String, String>> linked_props = ((Displayable) dbo).getLinkedProperties();
                        if (null != linked_props) {
                            for (final Map.Entry<Displayable, Map<String, String>> e : linked_props.entrySet()) {
                                for (final Map.Entry<String, String> ee : e.getValue().entrySet()) {
                                    if (pat.matcher(ee.getKey()).matches() || pat.matcher(ee.getValue()).matches()) {
                                        matched = true;
                                        txt = ee.getKey() + " => " + e.getValue() + " [linked property]";
                                        break;
                                    }
                                }
                            }
                        }
                    }
                }
                if (!matched && dbo instanceof Tree<?>) {
                    // search Node tags
                    Node<?> root = ((Tree<?>) dbo).getRoot();
                    if (null == root)
                        continue;
                    for (final Node<?> nd : root.getSubtreeNodes()) {
                        Set<Tag> tags = nd.getTags();
                        if (null == tags)
                            continue;
                        for (final Tag tag : tags) {
                            if (pat.matcher(tag.toString()).matches()) {
                                v_obs.add(dbo);
                                v_txt.add(new StringBuilder(tag.toString()).append(" (").append(null == meaningful_title ? dbo.toString() : meaningful_title).append(')').toString());
                                v_co.add(createCoordinate((Tree<?>) dbo, nd));
                            }
                        }
                    }
                    // all added if any
                    continue;
                }
                if (!matched)
                    continue;
                // txt = txt.length() > 30 ? txt.substring(0, 27) + "..." :
                // txt;
                v_obs.add(dbo);
                v_txt.add(txt);
                v_co.add(co);
            }
            if (0 == v_obs.size()) {
                Utils.showMessage("Nothing found.");
                return;
            }
            final JPanel result = new JPanel();
            GridBagLayout gb = new GridBagLayout();
            result.setLayout(gb);
            GridBagConstraints c = new GridBagConstraints();
            c.anchor = GridBagConstraints.NORTHWEST;
            c.fill = GridBagConstraints.HORIZONTAL;
            c.insets = new Insets(5, 10, 5, 10);
            String xml = "";
            if (project.getLoader() instanceof FSLoader) {
                String path = ((FSLoader) project.getLoader()).getProjectXMLPath();
                if (null != path) {
                    xml = " [" + new File(path).getName() + "]";
                }
            }
            JLabel projectTitle = new JLabel(project.getTitle() + xml);
            gb.setConstraints(projectTitle, c);
            result.add(projectTitle);
            c.insets = new Insets(0, 0, 0, 0);
            JPanel padding = new JPanel();
            c.weightx = 1;
            gb.setConstraints(padding, c);
            result.add(padding);
            c.gridy = 1;
            c.gridwidth = 2;
            c.fill = GridBagConstraints.BOTH;
            c.weighty = 1;
            JScrollPane jsp = makeTable(new DisplayableTableModel(v_obs, v_txt, v_co), project);
            gb.setConstraints(jsp, c);
            result.add(jsp);
            search_tabs.addTab(typed_pattern, result);
            search_tabs.setSelectedComponent(result);
            synchronized (tabMap) {
                List<JPanel> cs = tabMap.get(project);
                if (null == cs) {
                    cs = new ArrayList<JPanel>();
                    tabMap.put(project, cs);
                }
                cs.add(result);
            }
        }
    }, project);
}
Also used : JPanel(javax.swing.JPanel) GridBagConstraints(java.awt.GridBagConstraints) Set(java.util.Set) HashSet(java.util.HashSet) LayerSet(ini.trakem2.display.LayerSet) Insets(java.awt.Insets) GridBagLayout(java.awt.GridBagLayout) Node(ini.trakem2.display.Node) ArrayList(java.util.ArrayList) DBObject(ini.trakem2.persistence.DBObject) AreaTree(ini.trakem2.display.AreaTree) Tree(ini.trakem2.display.Tree) AreaList(ini.trakem2.display.AreaList) List(java.util.List) ArrayList(java.util.ArrayList) Vector(java.util.Vector) JScrollPane(javax.swing.JScrollPane) Pattern(java.util.regex.Pattern) Displayable(ini.trakem2.display.Displayable) ZDisplayable(ini.trakem2.display.ZDisplayable) JLabel(javax.swing.JLabel) Project(ini.trakem2.Project) FSLoader(ini.trakem2.persistence.FSLoader) Coordinate(ini.trakem2.display.Coordinate) Tag(ini.trakem2.display.Tag) Map(java.util.Map) HashMap(java.util.HashMap) File(java.io.File)

Example 20 with Coordinate

use of ini.trakem2.display.Coordinate in project TrakEM2 by trakem2.

the class AreaUtils method generateTriangles.

/**
 * Expects areas in local coordinates to the Displayable @param d.
 *  @param d
 *  @param scale The scaling of the entire universe, to limit the overall box
 *  @param resample_ The optimization parameter for marching cubes (i.e. a value of 2 will scale down to half, then apply marching cubes, then scale up by 2 the vertices coordinates).
 *  @param areas
 *  @return The List of triangles involved, specified as three consecutive vertices. A list of Point3f vertices.
 */
public static List<Point3f> generateTriangles(final Displayable d, final double scale, final int resample_, final Map<Layer, Area> areas) {
    // in the LayerSet, layers are ordered by Z already.
    try {
        int n = areas.size();
        if (0 == n)
            return null;
        final int resample;
        if (resample_ <= 0) {
            resample = 1;
            Utils.log2("Fixing zero or negative resampling value to 1.");
        } else
            resample = resample_;
        final LayerSet layer_set = d.getLayerSet();
        final AffineTransform aff = d.getAffineTransformCopy();
        final Rectangle r = d.getBoundingBox(null);
        // remove translation from a copy of the Displayable's AffineTransform
        final AffineTransform at_translate = new AffineTransform();
        at_translate.translate(-r.x, -r.y);
        aff.preConcatenate(at_translate);
        // incorporate resampling scaling into the transform
        final AffineTransform atK = new AffineTransform();
        // Utils.log("resample: " + resample + "  scale: " + scale);
        // 'scale' is there to limit gigantic universes
        final double K = (1.0 / resample) * scale;
        atK.scale(K, K);
        aff.preConcatenate(atK);
        final Calibration cal = layer_set.getCalibrationCopy();
        // Find first layer, compute depth, and fill in the depth vs area map
        Layer first_layer = null, last_layer = null;
        final int w = (int) Math.ceil(r.width * K);
        final int h = (int) Math.ceil(r.height * K);
        int depth = 0;
        final Map<Integer, Area> ma = new HashMap<Integer, Area>();
        for (final Layer la : layer_set.getLayers()) {
            // layers sorted by Z ASC
            final Area area = areas.get(la);
            if (null != area) {
                ma.put(depth, area);
                if (null == first_layer) {
                    first_layer = la;
                }
                // Utils.log("area at depth " + depth + " for layer " + la);
                depth++;
                n--;
            } else if (0 != depth) {
                // Utils.log("Empty area at depth " + depth);
                // an empty layer
                depth++;
            }
            if (0 == n) {
                last_layer = la;
                // no more areas to paint
                break;
            }
        }
        if (0 == depth) {
            Utils.log("ERROR could not find any areas for " + d);
            return null;
        }
        if (0 != n) {
            Utils.log("WARNING could not find all areas for " + d);
        }
        // No zero-padding: Marching Cubes now can handle edges
        final ShapeList<ByteType> shapeList = new ShapeListCached<ByteType>(new int[] { w, h, depth }, new ByteType(), 32);
        final Image<ByteType> shapeListImage = new Image<ByteType>(shapeList, shapeList.getBackground(), "ShapeListContainer");
        // 255 or -1 don't work !? So, giving the highest value (127) that is both a byte and an int.
        final ByteType intensity = new ByteType((byte) 127);
        for (final Map.Entry<Integer, Area> e : ma.entrySet()) {
            Area a = e.getValue();
            if (!aff.isIdentity()) {
                a = M.areaInIntsByRounding(a.createTransformedArea(aff));
            }
            shapeList.addShape(a, intensity, new int[] { e.getKey() });
        }
        // debug:
        // ImagePlus imp = ImageJFunctions.displayAsVirtualStack(shapeListImage);
        // imp.getProcessor().setMinAndMax( 0, 255 );
        // imp.show();
        // Utils.log2("Using imglib Shape List Image Container");
        // Now marching cubes
        // origins at 0,0,0: uncalibrated
        final List<Point3f> list = new MCTriangulator().getTriangles(shapeListImage, 1, new float[3]);
        // The list of triangles has coordinates:
        // - in x,y: in pixels, scaled by K = (1 / resample) * scale,
        // translated by r.x, r.y (the top-left coordinate of this AreaList bounding box)
        // - in z: in stack slice indices
        // So all x,y,z must be corrected in x,y and z of the proper layer
        // final double offset = first_layer.getZ();
        final int i_first_layer = layer_set.indexOf(first_layer);
        // The x,y translation to correct each point by:
        final float dx = (float) (r.x * scale * cal.pixelWidth);
        final float dy = (float) (r.y * scale * cal.pixelHeight);
        // Correct x,y by resampling and calibration, but not scale
        // scale is already in the pixel coordinates
        final float rsw = (float) (resample * cal.pixelWidth);
        final float rsh = (float) (resample * cal.pixelHeight);
        // no resampling in Z. and Uses pixelWidth, not pixelDepth.
        final double sz = scale * cal.pixelWidth;
        // debug:
        /*
			// which p.z types exist?
			final TreeSet<Float> ts = new TreeSet<Float>();
			for (final Iterator it = list.iterator(); it.hasNext(); ) {
				ts.add(((Point3f)it.next()).z);
			}
			for (final Float pz : ts) Utils.log2("A z: " + pz);
			*/
        // debug: How many different Z?
        /*
			HashSet<Float> zs = new HashSet<Float>();
			for (Point3f p : list) {
				zs.add(p.z);
			}
			ArrayList<Float> a = new ArrayList<Float>(zs);
			java.util.Collections.sort(a);
			for (Float f : a) {
				Utils.log("f: " + f);
			}
			*/
        // Utils.log2("Number of slices: " + imp.getNSlices());
        // Fix all points:
        // Read from list, modify and put into verts
        // and don't modify it if the verts already has it (it's just coincident)
        final Point3f[] verts = new Point3f[list.size()];
        // Utils.log("number of verts: " + verts.length + " mod 3: " + (verts.length % 3));
        final TreeMap<Integer, Point3f> output = new TreeMap<Integer, Point3f>();
        // The first section generates vertices at -1 and 0
        // The last section generates them at last_section_index and last_section_index +1
        // Capture from -1 to 0
        fix3DPoints(list, output, verts, first_layer.getZ(), 0, -1, dx, dy, rsw, rsh, sz, 1);
        int slice_index = 0;
        for (final Layer la : layer_set.getLayers().subList(i_first_layer, i_first_layer + depth)) {
            // If layer is empty, continue
            /* // YEAH don't! At least the immediate next layer would have points, like the extra Z level after last layer, to account for the thickness of the layer!
				if (empty_layers.contains(la)) {
					slice_index++;
					continue;
				}
				*/
            fix3DPoints(list, output, verts, la.getZ(), la.getThickness(), slice_index, dx, dy, rsw, rsh, sz, 1);
            slice_index++;
        }
        // Do the last layer again. The last layer has two Z planes in which it has pixels:
        try {
            // Capture from last_section_index to last_section_index+1, inclusive
            fix3DPoints(list, output, verts, last_layer.getZ() + last_layer.getThickness(), 0, slice_index, dx, dy, rsw, rsh, sz, 2);
        } catch (final Exception ee) {
            IJError.print(ee);
        }
        // Handle potential errors:
        if (0 != list.size() - output.size()) {
            Utils.log2("Unprocessed/unused points: " + (list.size() - output.size()));
            for (int i = 0; i < verts.length; i++) {
                if (null == verts[i]) {
                    final Point3f p = (Point3f) list.get(i);
                    Utils.log2("verts[" + i + "] = " + p.x + ", " + p.y + ", " + p.z + "  p.z as int: " + ((int) (p.z + 0.05f)));
                }
            }
            return new ArrayList<Point3f>(output.values());
        } else {
            return java.util.Arrays.asList(verts);
        }
    } catch (final Exception e) {
        e.printStackTrace();
    }
    return null;
}
Also used : HashMap(java.util.HashMap) Rectangle(java.awt.Rectangle) ArrayList(java.util.ArrayList) ShapeListCached(mpicbg.imglib.container.shapelist.ShapeListCached) ByteType(mpicbg.imglib.type.numeric.integer.ByteType) Image(mpicbg.imglib.image.Image) Point3f(org.scijava.vecmath.Point3f) LayerSet(ini.trakem2.display.LayerSet) Calibration(ij.measure.Calibration) TreeMap(java.util.TreeMap) Layer(ini.trakem2.display.Layer) ExecutionException(java.util.concurrent.ExecutionException) Area(java.awt.geom.Area) AffineTransform(java.awt.geom.AffineTransform) HashMap(java.util.HashMap) Map(java.util.Map) TreeMap(java.util.TreeMap)

Aggregations

ArrayList (java.util.ArrayList)14 Layer (ini.trakem2.display.Layer)10 HashMap (java.util.HashMap)10 Rectangle (java.awt.Rectangle)9 Patch (ini.trakem2.display.Patch)8 AffineTransform (java.awt.geom.AffineTransform)7 Map (java.util.Map)7 Displayable (ini.trakem2.display.Displayable)6 NoninvertibleTransformException (java.awt.geom.NoninvertibleTransformException)6 TreeMap (java.util.TreeMap)6 Worker (ini.trakem2.utils.Worker)5 Area (java.awt.geom.Area)5 HashSet (java.util.HashSet)5 GenericDialog (ij.gui.GenericDialog)4 Calibration (ij.measure.Calibration)4 ByteProcessor (ij.process.ByteProcessor)4 File (java.io.File)4 Future (java.util.concurrent.Future)4 ImagePlus (ij.ImagePlus)3 ColorProcessor (ij.process.ColorProcessor)3