Examples with Watershed3D mcib3d.image3d.regionGrowing.Watershed3D used on opensource projects

Example 1 with Watershed3D

use of mcib3d.image3d.regionGrowing.Watershed3D in project mcib3d-core by mcib3d.

the class Segment3DSpots method splitSpotWatershed.

public static Object3DVoxels[] splitSpotWatershed(Object3D obj, float rad, float dist) {
    ImageInt seg = obj.createSegImage(0, 0, 0, obj.getXmax(), obj.getYmax(), obj.getZmax(), 255);
    // seg.show();
    ImagePlus segplus = seg.getImagePlus();
    segplus.setCalibration(obj.getCalibration());
    // return
    Object3DVoxels[] res = null;
    try {
        int cpus = ThreadUtil.getNbCpus();
        // FIXME variable multithread
        ImageFloat edt3d = EDT.run(seg, 1f, false, cpus);
        // 3D filtering of the edt to remove small local maxima
        ImageFloat edt3dseeds = FastFilters3D.filterFloatImage(edt3d, FastFilters3D.MEAN, 1, 1, 1, cpus, false);
        // edt3d.showDuplicate("edt");
        // ImageStack localMax = FastFilters3D.filterFloatImageStack(edt3d.getImageStack(), FastFilters3D.MAXLOCAL, rad, rad, rad, cpus, false);
        ImageFloat maxlocal3d = FastFilters3D.filterFloatImage(edt3dseeds, FastFilters3D.MAXLOCAL, rad, rad, rad, cpus, false);
        // maxlocal3d.show("max local");
        ArrayList<Voxel3D> locals = obj.listVoxels(maxlocal3d, 0);
        int nb = locals.size();
        // IJ.log("nb=" + nb);
        if (nb < 2) {
            return null;
        }
        // look for most far apart local maxima
        int i1 = 0;
        int i2 = 0;
        double dmax = 0.0;
        double d;
        for (int i = 0; i < nb; i++) {
            for (int j = i + 1; j < nb; j++) {
                d = locals.get(i).distance(locals.get(j));
                if (d > dmax) {
                    dmax = d;
                    i1 = i;
                    i2 = j;
                }
            }
        }
        // perform a 2-means clustering
        double cx1 = 0;
        double cy1 = 0;
        double cz1 = 0;
        double cx2 = 0;
        double cy2 = 0;
        double cz2 = 0;
        double d1;
        double d2;
        Voxel3D PP1 = new Voxel3D(locals.get(i1).getX(), locals.get(i1).getY(), locals.get(i1).getZ(), 1);
        Voxel3D PP2 = new Voxel3D(locals.get(i2).getX(), locals.get(i2).getY(), locals.get(i2).getZ(), 2);
        Voxel3D P1 = new Voxel3D(cx1, cy1, cz1, 0);
        Voxel3D P2 = new Voxel3D(cx2, cy2, cz2, 0);
        int nb1, nb2;
        int nbite = 0;
        while ((nb > 2) && ((P1.distance(PP1) > 1) || (P2.distance(PP2) > 1)) && (nbite < 100)) {
            nbite++;
            cx1 = 0;
            cy1 = 0;
            cx2 = 0;
            cy2 = 0;
            cz1 = 0;
            cz2 = 0;
            nb1 = 0;
            nb2 = 0;
            P1.setX(PP1.getX());
            P1.setY(PP1.getY());
            P1.setZ(PP1.getZ());
            P2.setX(PP2.getX());
            P2.setY(PP2.getY());
            P2.setZ(PP2.getZ());
            for (Voxel3D local : locals) {
                d1 = P1.distance(local);
                d2 = P2.distance(local);
                if (d1 < d2) {
                    cx1 += local.getX();
                    cy1 += local.getY();
                    cz1 += local.getZ();
                    nb1++;
                } else {
                    cx2 += local.getX();
                    cy2 += local.getY();
                    cz2 += local.getZ();
                    nb2++;
                }
            }
            cx1 /= nb1;
            cy1 /= nb1;
            cx2 /= nb2;
            cy2 /= nb2;
            cz1 /= nb1;
            cz2 /= nb2;
            PP1.setX(cx1);
            PP1.setY(cy1);
            PP1.setZ(cz1);
            PP2.setX(cx2);
            PP2.setY(cy2);
            PP2.setZ(cz2);
        }
        // check minimal distances
        double distPP = PP1.distance(PP2);
        IJ.log("Centers found for split PP1=" + PP1 + " PP2=" + PP2 + " distance " + distPP);
        if (distPP < dist) {
            return null;
        }
        // find closest max local to two barycenters, in case barycenters outside object
        Voxel3D PP1closest = new Voxel3D(locals.get(0));
        Voxel3D PP2closest = new Voxel3D(locals.get(0));
        double dist1 = PP1.distanceSquare(PP1closest);
        double dist2 = PP2.distanceSquare(PP2closest);
        for (Voxel3D local : locals) {
            double dd1 = PP1.distanceSquare(local);
            double dd2 = PP2.distanceSquare(local);
            if (dd1 < dist1) {
                dist1 = dd1;
                PP1closest = local;
            }
            if (dd2 < dist2) {
                dist2 = dd2;
                PP2closest = local;
            }
        }
        ImageInt seeds = new ImageShort("seeds", seg.sizeX, seg.sizeY, seg.sizeZ);
        seeds.setPixel(PP1closest.getRoundX(), PP1closest.getRoundY(), PP1closest.getRoundZ(), 255);
        seeds.setPixel(PP2closest.getRoundX(), PP2closest.getRoundY(), PP2closest.getRoundZ(), 255);
        // seeds.show();
        // Watershed3D_old wat = new Watershed3D_old(edt3d, seeds, 0, 0);
        // ImageInt wat2 = wat.getWatershedImage3D();
        // ImageHandler edt16 = edt3d.convertToShort(true);
        Watershed3D wat = new Watershed3D(edt3d, seeds, 0, 0);
        ImageInt wat2 = wat.getWatershedImage3D();
        // wat2.show();
        // in watershed label starts at 1
        Object3DVoxels ob1 = new Object3DVoxels(wat2, 1);
        Object3DVoxels ob2 = new Object3DVoxels(wat2, 2);
        // IJ.log("split1="+ob1+" split2="+ob2);
        // translate objects if needed by miniseg
        // ob1.translate(seg.offsetX, seg.offsetY, seg.offsetZ);
        // new ImagePlus("wat", wat2.getStack()).show();
        res = new Object3DVoxels[2];
        res[0] = ob1;
        res[1] = ob2;
    } catch (Exception e) {
        IJ.log("Exception EDT " + e);
    }
    return res;
}

Example 2 with Watershed3D

use of mcib3d.image3d.regionGrowing.Watershed3D in project mcib3d-core by mcib3d.

the class Segment3DSpots method computeWatershed.

private void computeWatershed() {
    // watershed is used to separate spots, not for segmentation
    // so based on edt of background
    ImageByte seedsTh = seedsImage.thresholdAboveExclusive(seedsThreshold);
    ImageFloat edt = EDT.run(seedsTh, 0, (float) rawImage.getScaleXY(), (float) rawImage.getScaleZ(), true, 0);
    ImageShort edt16 = edt.convertToShort(true);
    edt16.invert();
    Watershed3D wat3D = new Watershed3D(edt16, seedsImage, 0, 0);
    watershedImage = wat3D.getWatershedImage3D();
}