Examples with ExtendedRandomAccessibleInterval net.imglib2.view.ExtendedRandomAccessibleInterval used on opensource projects

Example 1 with ExtendedRandomAccessibleInterval

use of net.imglib2.view.ExtendedRandomAccessibleInterval in project imagej-ops by imagej.

the class LocalThresholdIntegral method getIntegralImage.

/**
 * Computes integral images of a given order and extends them such that
 * {@link IntegralMean} et al work with them.
 *
 * @param input The RAI for which an integral image is computed
 * @param order
 * @return An extended integral image for the input RAI
 */
@SuppressWarnings({ "unchecked", "rawtypes" })
private RandomAccessibleInterval<RealType> getIntegralImage(final RandomAccessibleInterval<I> input, final int order) {
    ExtendedRandomAccessibleInterval<I, RandomAccessibleInterval<I>> extendedInput = Views.extend(input, outOfBoundsFactory);
    FinalInterval expandedInterval = Intervals.expand(input, shape.getSpan() - 1);
    IntervalView<I> offsetInterval2 = Views.offsetInterval(extendedInput, expandedInterval);
    RandomAccessibleInterval<RealType> img = null;
    switch(order) {
        case 1:
            img = (RandomAccessibleInterval) integralImgOp.calculate(offsetInterval2);
            break;
        case 2:
            img = (RandomAccessibleInterval) squareIntegralImgOp.calculate(offsetInterval2);
            break;
    }
    img = addLeadingZeros(img);
    return img;
}

Example 2 with ExtendedRandomAccessibleInterval

use of net.imglib2.view.ExtendedRandomAccessibleInterval in project imagej-ops by imagej.

the class DefaultMarchingCubes method calculate.

@SuppressWarnings({ "unchecked" })
@Override
public DefaultMesh calculate(final RandomAccessibleInterval<T> input) {
    DefaultMesh output = new DefaultMesh();
    ExtendedRandomAccessibleInterval<T, RandomAccessibleInterval<T>> extended = Views.extendValue(input, (T) new BoolType(false));
    Cursor<T> c = Views.interval(extended, new FinalInterval(new long[] { input.min(0) - 1, input.min(1) - 1, input.min(2) - 1 }, new long[] { input.max(0) + 1, input.max(1) + 1, input.max(2) + 1 })).localizingCursor();
    while (c.hasNext()) {
        c.next();
        int cursorX = c.getIntPosition(0);
        int cursorY = c.getIntPosition(1);
        int cursorZ = c.getIntPosition(2);
        Cursor<T> cu = getCube(extended, cursorX, cursorY, cursorZ);
        int i = 0;
        double[] vertex_values = new double[8];
        while (cu.hasNext()) {
            vertex_values[i++] = (cu.next().get()) ? 1 : 0;
        }
        // 6------7
        // /| /|
        // 2-----3 |
        // | 4---|-5
        // |/ |/
        // 0-----1
        vertex_values = mapFlatIterableToLookUpCube(vertex_values);
        // 4------5
        // /| /|
        // 7-----6 |
        // | 0---|-1
        // |/ |/
        // 3-----2
        int cubeindex = getCubeIndex(vertex_values);
        if (EDGE_TABLE[cubeindex] != 0) {
            int[] p0 = new int[] { 0 + cursorX, 0 + cursorY, 1 + cursorZ };
            int[] p1 = new int[] { 1 + cursorX, 0 + cursorY, 1 + cursorZ };
            int[] p2 = new int[] { 1 + cursorX, 0 + cursorY, 0 + cursorZ };
            int[] p3 = new int[] { 0 + cursorX, 0 + cursorY, 0 + cursorZ };
            int[] p4 = new int[] { 0 + cursorX, 1 + cursorY, 1 + cursorZ };
            int[] p5 = new int[] { 1 + cursorX, 1 + cursorY, 1 + cursorZ };
            int[] p6 = new int[] { 1 + cursorX, 1 + cursorY, 0 + cursorZ };
            int[] p7 = new int[] { 0 + cursorX, 1 + cursorY, 0 + cursorZ };
            double[][] vertlist = new double[12][];
            /* Find the vertices where the surface intersects the cube */
            if (0 != (EDGE_TABLE[cubeindex] & 1)) {
                vertlist[0] = interpolatePoint(p0, p1, vertex_values[0], vertex_values[1]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 2)) {
                vertlist[1] = interpolatePoint(p1, p2, vertex_values[1], vertex_values[2]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 4)) {
                vertlist[2] = interpolatePoint(p2, p3, vertex_values[2], vertex_values[3]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 8)) {
                vertlist[3] = interpolatePoint(p3, p0, vertex_values[3], vertex_values[0]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 16)) {
                vertlist[4] = interpolatePoint(p4, p5, vertex_values[4], vertex_values[5]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 32)) {
                vertlist[5] = interpolatePoint(p5, p6, vertex_values[5], vertex_values[6]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 64)) {
                vertlist[6] = interpolatePoint(p6, p7, vertex_values[6], vertex_values[7]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 128)) {
                vertlist[7] = interpolatePoint(p7, p4, vertex_values[7], vertex_values[4]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 256)) {
                vertlist[8] = interpolatePoint(p0, p4, vertex_values[0], vertex_values[4]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 512)) {
                vertlist[9] = interpolatePoint(p1, p5, vertex_values[1], vertex_values[5]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 1024)) {
                vertlist[10] = interpolatePoint(p2, p6, vertex_values[2], vertex_values[6]);
            }
            if (0 != (EDGE_TABLE[cubeindex] & 2048)) {
                vertlist[11] = interpolatePoint(p3, p7, vertex_values[3], vertex_values[7]);
            }
            /* Create the triangle */
            for (i = 0; TRIANGLE_TABLE[cubeindex][i] != -1; i += 3) {
                TriangularFacet face = new TriangularFacet(new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][0], vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][1], vertlist[TRIANGLE_TABLE[cubeindex][i + 2]][2]), new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][0], vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][1], vertlist[TRIANGLE_TABLE[cubeindex][i + 1]][2]), new Vertex(vertlist[TRIANGLE_TABLE[cubeindex][i]][0], vertlist[TRIANGLE_TABLE[cubeindex][i]][1], vertlist[TRIANGLE_TABLE[cubeindex][i]][2]));
                face.getArea();
                output.addFace(face);
            }
        }
    }
    return output;
}

Example 3 with ExtendedRandomAccessibleInterval

use of net.imglib2.view.ExtendedRandomAccessibleInterval in project imagej-ops by imagej.

the class Watershed method compute.

@Override
public void compute(final RandomAccessibleInterval<T> in, final ImgLabeling<Integer, IntType> out) {
    final RandomAccess<T> raIn = in.randomAccess();
    RandomAccess<B> raMask = null;
    if (mask != null) {
        raMask = mask.randomAccess();
    }
    // stores the size of each dimension
    final long[] dimensSizes = new long[in.numDimensions()];
    in.dimensions(dimensSizes);
    // calculates the number of points in the n-d space
    long numPixels = Intervals.numElements(in);
    // the pixels indices are stored in an array, which is sorted depending
    // on the pixel values
    final List<Long> imiList = new ArrayList<>();
    if (mask != null) {
        final Cursor<Void> c = Regions.iterable(mask).localizingCursor();
        while (c.hasNext()) {
            c.next();
            imiList.add(IntervalIndexer.positionToIndex(c, in));
        }
    } else {
        for (long i = 0; i < numPixels; i++) {
            imiList.add(i);
        }
    }
    final Long[] imi = imiList.toArray(new Long[imiList.size()]);
    /*
		 * Sort the pixels of imi in the increasing order of their grey value
		 * (only the pixel indices are stored)
		 */
    Arrays.sort(imi, new Comparator<Long>() {

        @Override
        public int compare(final Long o1, final Long o2) {
            IntervalIndexer.indexToPosition(o1, in, raIn);
            final T value = raIn.get().copy();
            IntervalIndexer.indexToPosition(o2, in, raIn);
            return value.compareTo(raIn.get());
        }
    });
    // lab and dist store the values calculated after each phase
    final RandomAccessibleInterval<IntType> lab = ops().create().img(in, new IntType());
    // extend border to be able to do a quick check, if a voxel is inside
    final ExtendedRandomAccessibleInterval<IntType, RandomAccessibleInterval<IntType>> labExt = Views.extendBorder(lab);
    final OutOfBounds<IntType> raLab = labExt.randomAccess();
    final RandomAccessibleInterval<IntType> dist = ops().create().img(in, new IntType());
    final RandomAccess<IntType> raDist = dist.randomAccess();
    // initial values
    for (final IntType pixel : Views.flatIterable(lab)) {
        pixel.set(INIT);
    }
    int current_label = 0;
    int current_dist;
    final ArrayList<Long> fifo = new ArrayList<>();
    // RandomAccess for Neighborhoods
    final Shape shape;
    if (useEightConnectivity) {
        shape = new RectangleShape(1, true);
    } else {
        shape = new DiamondShape(1);
    }
    final RandomAccessible<Neighborhood<T>> neighborhoods = shape.neighborhoodsRandomAccessible(in);
    final RandomAccess<Neighborhood<T>> raNeighbor = neighborhoods.randomAccess();
    /*
		 * Start flooding
		 */
    for (int j = 0; j < imi.length; j++) {
        IntervalIndexer.indexToPosition(imi[j], in, raIn);
        final T actualH = raIn.get().copy();
        int i = j;
        while (actualH.compareTo(raIn.get()) == 0) {
            final long p = imi[i];
            IntervalIndexer.indexToPosition(p, in, raIn);
            raLab.setPosition(raIn);
            raLab.get().set(MASK);
            raNeighbor.setPosition(raIn);
            final Cursor<T> neighborHood = raNeighbor.get().cursor();
            while (neighborHood.hasNext()) {
                neighborHood.fwd();
                raLab.setPosition(neighborHood);
                if (!raLab.isOutOfBounds()) {
                    final int f = raLab.get().get();
                    if ((f > 0) || (f == WSHED)) {
                        raDist.setPosition(raIn);
                        raDist.get().set(1);
                        fifo.add(p);
                        break;
                    }
                }
            }
            i++;
            if (i == imi.length) {
                break;
            }
            IntervalIndexer.indexToPosition(imi[i], in, raIn);
        }
        current_dist = 1;
        // add fictitious pixel
        fifo.add(-1l);
        while (true) {
            long p = fifo.remove(0);
            if (p == -1) {
                if (fifo.isEmpty()) {
                    break;
                }
                fifo.add(-1l);
                current_dist++;
                p = fifo.remove(0);
            }
            IntervalIndexer.indexToPosition(p, in, raNeighbor);
            final Cursor<T> neighborHood = raNeighbor.get().cursor();
            raLab.setPosition(raNeighbor);
            int labp = raLab.get().get();
            final long[] posNeighbor = new long[neighborHood.numDimensions()];
            while (neighborHood.hasNext()) {
                neighborHood.fwd();
                neighborHood.localize(posNeighbor);
                raLab.setPosition(posNeighbor);
                if (!raLab.isOutOfBounds()) {
                    raDist.setPosition(posNeighbor);
                    final int labq = raLab.get().get();
                    final int distq = raDist.get().get();
                    if ((distq < current_dist) && ((labq > 0) || (labq == WSHED))) {
                        // the watersheds
                        if (labq > 0) {
                            if ((labp == MASK) || (labp == WSHED)) {
                                labp = labq;
                            } else {
                                if (labp != labq) {
                                    labp = WSHED;
                                }
                            }
                        } else {
                            if (labp == MASK) {
                                labp = WSHED;
                            }
                        }
                        raLab.setPosition(raNeighbor);
                        raLab.get().set(labp);
                    } else {
                        if ((labq == MASK) && (distq == 0)) {
                            raDist.setPosition(posNeighbor);
                            raDist.get().set(current_dist + 1);
                            fifo.add(IntervalIndexer.positionToIndex(posNeighbor, dimensSizes));
                        }
                    }
                }
            }
        }
        // checks if new minima have been discovered
        IntervalIndexer.indexToPosition(imi[j], in, raIn);
        i = j;
        while (actualH.compareTo(raIn.get()) == 0) {
            final long p = imi[i];
            IntervalIndexer.indexToPosition(p, dist, raDist);
            // the distance associated with p is reseted to 0
            raDist.get().set(0);
            raLab.setPosition(raDist);
            if (raLab.get().get() == MASK) {
                current_label++;
                fifo.add(p);
                raLab.get().set(current_label);
                while (!fifo.isEmpty()) {
                    final long q = fifo.remove(0);
                    IntervalIndexer.indexToPosition(q, in, raNeighbor);
                    final Cursor<T> neighborHood = raNeighbor.get().cursor();
                    final long[] posNeighbor = new long[neighborHood.numDimensions()];
                    while (neighborHood.hasNext()) {
                        neighborHood.fwd();
                        neighborHood.localize(posNeighbor);
                        raLab.setPosition(posNeighbor);
                        if (!raLab.isOutOfBounds()) {
                            final long r = IntervalIndexer.positionToIndex(posNeighbor, dimensSizes);
                            if (raLab.get().get() == MASK) {
                                fifo.add(r);
                                raLab.get().set(current_label);
                            }
                        }
                    }
                }
            }
            i++;
            if (i == imi.length) {
                break;
            }
            IntervalIndexer.indexToPosition(imi[i], in, raIn);
        }
        j = i - 1;
    }
    /*
		 * Draw output and remove as the case may be the watersheds
		 */
    final Cursor<LabelingType<Integer>> cursorOut = out.cursor();
    while (cursorOut.hasNext()) {
        cursorOut.fwd();
        boolean maskValue = true;
        if (mask != null) {
            raMask.setPosition(cursorOut);
            if (!raMask.get().get()) {
                maskValue = false;
            }
        }
        raLab.setPosition(cursorOut);
        if (!maskValue) {
            cursorOut.get().clear();
        } else {
            if (!drawWatersheds && raLab.get().get() == WSHED) {
                raNeighbor.setPosition(cursorOut);
                final Cursor<T> neighborHood = raNeighbor.get().cursor();
                int newLab = WSHED;
                while (neighborHood.hasNext()) {
                    neighborHood.fwd();
                    raLab.setPosition(neighborHood);
                    if (!raLab.isOutOfBounds()) {
                        newLab = raLab.get().get();
                        if (newLab > WSHED) {
                            break;
                        }
                    }
                }
                if (newLab == WSHED) {
                    cursorOut.get().clear();
                } else {
                    cursorOut.get().add(newLab);
                }
            } else {
                cursorOut.get().add(raLab.get().get());
            }
        }
    }
    /*
		 * Merge already present labels before calculation of watershed
		 */
    if (out() != null) {
        final Cursor<LabelingType<Integer>> cursor = out().cursor();
        final RandomAccess<LabelingType<Integer>> raOut = out.randomAccess();
        while (cursor.hasNext()) {
            cursor.fwd();
            raOut.setPosition(cursor);
            final List<Integer> labels = new ArrayList<>();
            cursor.get().iterator().forEachRemaining(labels::add);
            raOut.get().addAll(labels);
        }
    }
}

Example 4 with ExtendedRandomAccessibleInterval

use of net.imglib2.view.ExtendedRandomAccessibleInterval in project imagej-ops by imagej.

the class WatershedSeeded method compute.

@SuppressWarnings("unchecked")
@Override
public void compute(final RandomAccessibleInterval<T> in, final ImgLabeling<Integer, IntType> out) {
    // extend border to be able to do a quick check, if a voxel is inside
    final LabelingType<Integer> oustide = out.firstElement().copy();
    oustide.clear();
    oustide.add(OUTSIDE);
    final ExtendedRandomAccessibleInterval<LabelingType<Integer>, ImgLabeling<Integer, IntType>> outExt = Views.extendValue(out, oustide);
    final OutOfBounds<LabelingType<Integer>> raOut = outExt.randomAccess();
    // if no mask provided, set the mask to the whole image
    if (mask == null) {
        mask = (RandomAccessibleInterval<B>) ops().create().img(in, new BitType());
        for (B b : Views.flatIterable(mask)) {
            b.set(true);
        }
    }
    // initialize output labels
    final Cursor<B> maskCursor = Views.flatIterable(mask).cursor();
    while (maskCursor.hasNext()) {
        maskCursor.fwd();
        if (maskCursor.get().get()) {
            raOut.setPosition(maskCursor);
            raOut.get().clear();
            raOut.get().add(INIT);
        }
    }
    // RandomAccess for Mask, Seeds and Neighborhoods
    final RandomAccess<B> raMask = mask.randomAccess();
    final RandomAccess<LabelingType<Integer>> raSeeds = seeds.randomAccess();
    final Shape shape;
    if (useEightConnectivity) {
        shape = new RectangleShape(1, true);
    } else {
        shape = new DiamondShape(1);
    }
    final RandomAccessible<Neighborhood<T>> neighborhoods = shape.neighborhoodsRandomAccessible(in);
    final RandomAccess<Neighborhood<T>> raNeigh = neighborhoods.randomAccess();
    /*
		 * Carry over the seeding points to the new label and adds them to a
		 * voxel priority queue
		 */
    final PriorityQueue<WatershedVoxel> pq = new PriorityQueue<>();
    // Only iterate seeds that are not excluded by the mask
    final IterableRegion<B> maskRegions = Regions.iterable(mask);
    final IterableInterval<LabelingType<Integer>> seedsMasked = Regions.sample(maskRegions, seeds);
    final Cursor<LabelingType<Integer>> cursorSeeds = seedsMasked.localizingCursor();
    while (cursorSeeds.hasNext()) {
        final Set<Integer> l = cursorSeeds.next();
        if (l.isEmpty()) {
            continue;
        }
        if (l.size() > 1) {
            throw new IllegalArgumentException("Seeds must have exactly one label!");
        }
        final Integer label = l.iterator().next();
        if (label < 0) {
            throw new IllegalArgumentException("Seeds must have positive integers as labels!");
        }
        raNeigh.setPosition(cursorSeeds);
        final Cursor<T> neighborhood = raNeigh.get().cursor();
        // Add unlabeled neighbors to priority queue
        while (neighborhood.hasNext()) {
            neighborhood.fwd();
            raSeeds.setPosition(neighborhood);
            raMask.setPosition(neighborhood);
            raOut.setPosition(neighborhood);
            final Integer labelNeigh = raOut.get().iterator().next();
            if (labelNeigh != INQUEUE && labelNeigh != OUTSIDE && !raOut.isOutOfBounds() && raMask.get().get() && raSeeds.get().isEmpty()) {
                raOut.setPosition(neighborhood);
                pq.add(new WatershedVoxel(IntervalIndexer.positionToIndex(neighborhood, in), neighborhood.get().getRealDouble()));
                raOut.get().clear();
                raOut.get().add(INQUEUE);
            }
        }
        // Overwrite label in output with the seed label
        raOut.setPosition(cursorSeeds);
        raOut.get().clear();
        raOut.get().add(label);
    }
    /*
		 * Pop the head of the priority queue, label and push all unlabeled
		 * neighbored pixels.
		 */
    // list to store neighbor labels
    final ArrayList<Integer> neighborLabels = new ArrayList<>();
    // list to store neighbor voxels
    final ArrayList<WatershedVoxel> neighborVoxels = new ArrayList<>();
    // iterate the queue
    final Point pos = new Point(in.numDimensions());
    while (!pq.isEmpty()) {
        IntervalIndexer.indexToPosition(pq.poll().getPos(), out, pos);
        // reset list of neighbor labels
        neighborLabels.clear();
        // reset list of neighbor voxels
        neighborVoxels.clear();
        // iterate the neighborhood of the pixel
        raNeigh.setPosition(pos);
        final Cursor<T> neighborhood = raNeigh.get().cursor();
        while (neighborhood.hasNext()) {
            neighborhood.fwd();
            // Unlabeled neighbors go into the queue if they are not there
            // yet
            raOut.setPosition(neighborhood);
            raMask.setPosition(raOut);
            if (!raOut.get().isEmpty()) {
                final Integer label = raOut.get().iterator().next();
                if (label == INIT && raMask.get().get()) {
                    neighborVoxels.add(new WatershedVoxel(IntervalIndexer.positionToIndex(neighborhood, out), neighborhood.get().getRealDouble()));
                } else {
                    if (label > WSHED && (!drawWatersheds || !neighborLabels.contains(label))) {
                        // store labels of neighbors in a list
                        neighborLabels.add(label);
                    }
                }
            }
        }
        if (drawWatersheds) {
            // if the neighbors of the extracted voxel that have already
            // been labeled
            // all have the same label, then the voxel is labeled with their
            // label.
            raOut.setPosition(pos);
            raOut.get().clear();
            if (neighborLabels.size() == 1) {
                raOut.get().add(neighborLabels.get(0));
                // list
                for (final WatershedVoxel v : neighborVoxels) {
                    IntervalIndexer.indexToPosition(v.getPos(), out, raOut);
                    raOut.get().clear();
                    raOut.get().add(INQUEUE);
                    pq.add(v);
                }
            } else if (neighborLabels.size() > 1)
                raOut.get().add(WSHED);
        } else {
            if (neighborLabels.size() > 0) {
                raOut.setPosition(pos);
                raOut.get().clear();
                // take the label which most of the neighbors have
                if (neighborLabels.size() > 2) {
                    final Map<Integer, Long> countLabels = neighborLabels.stream().collect(Collectors.groupingBy(e -> e, Collectors.counting()));
                    final Integer keyMax = Collections.max(countLabels.entrySet(), Comparator.comparingLong(Map.Entry::getValue)).getKey();
                    raOut.get().add(keyMax);
                } else {
                    raOut.get().add(neighborLabels.get(0));
                }
                // list
                for (final WatershedVoxel v : neighborVoxels) {
                    IntervalIndexer.indexToPosition(v.getPos(), out, raOut);
                    raOut.get().clear();
                    raOut.get().add(INQUEUE);
                    pq.add(v);
                }
            }
        }
    }
    /*
		 * Merge already present labels before calculation of watershed
		 */
    if (out() != null) {
        final Cursor<LabelingType<Integer>> cursor = out().cursor();
        while (cursor.hasNext()) {
            cursor.fwd();
            raOut.setPosition(cursor);
            final List<Integer> labels = new ArrayList<>();
            cursor.get().iterator().forEachRemaining(labels::add);
            raOut.get().addAll(labels);
        }
    }
}

Example 5 with ExtendedRandomAccessibleInterval

use of net.imglib2.view.ExtendedRandomAccessibleInterval in project imagej-ops by imagej.

the class Outline method neighbourhoodInterval.

/**
 * Creates a view that spans from (x-1, y-1, ... i-1) to (x+1, y+1, ... i+1)
 * around the given coordinates
 *
 * @param interval the space of the coordinates
 * @param coordinates coordinates (x, y, ... i)
 * @return a view of a neighbourhood in the space
 */
private IntervalView<B> neighbourhoodInterval(final ExtendedRandomAccessibleInterval<B, RandomAccessibleInterval<B>> interval, final long[] coordinates) {
    final int dimensions = interval.numDimensions();
    final BoundingBox box = new BoundingBox(dimensions);
    final long[] minBounds = Arrays.stream(coordinates).map(c -> c - 1).toArray();
    final long[] maxBounds = Arrays.stream(coordinates).map(c -> c + 1).toArray();
    box.update(minBounds);
    box.update(maxBounds);
    return Views.offsetInterval(interval, box);
}