Examples with Cursor net.imglib2.Cursor used on opensource projects

Example 61 with Cursor

use of net.imglib2.Cursor 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 62 with Cursor

use of net.imglib2.Cursor in project imagej-ops by imagej.

the class IntegralMean method compute.

@Override
public void compute(final RectangleNeighborhood<Composite<I>> input, final DoubleType output) {
    // computation according to
    // https://en.wikipedia.org/wiki/Summed_area_table
    final IntegralCursor<Composite<I>> cursor = new IntegralCursor<>(input);
    final int dimensions = input.numDimensions();
    // Compute \sum (-1)^{dim - ||cornerVector||_{1}} * I(x^{cornerVector})
    final DoubleType sum = new DoubleType();
    sum.setZero();
    // Convert from input to return type
    final Converter<I, DoubleType> conv = new RealDoubleConverter<>();
    final DoubleType valueAsDoubleType = new DoubleType();
    while (cursor.hasNext()) {
        final I value = cursor.next().get(0).copy();
        conv.convert(value, valueAsDoubleType);
        // Obtain the cursor position encoded as corner vector
        final int cornerInteger = cursor.getCornerRepresentation();
        // Determine if the value has to be added (factor==1) or subtracted
        // (factor==-1)
        final DoubleType factor = new DoubleType(Math.pow(-1.0d, dimensions - IntegralMean.norm(cornerInteger)));
        valueAsDoubleType.mul(factor);
        sum.add(valueAsDoubleType);
    }
    final int area = (int) Intervals.numElements(Intervals.expand(input, -1l));
    // Compute mean by dividing the sum divided by the number of elements
    // NB: Reuse DoubleType
    valueAsDoubleType.set(area);
    sum.div(valueAsDoubleType);
    output.set(sum);
}

Example 63 with Cursor

use of net.imglib2.Cursor in project imagej-ops by imagej.

the class IntegralSum method compute.

@Override
public void compute(final RectangleNeighborhood<I> input, final DoubleType output) {
    // computation according to
    // https://en.wikipedia.org/wiki/Summed_area_table
    final IntegralCursor<I> cursor = new IntegralCursor<>(input);
    final int dimensions = input.numDimensions();
    // Compute \sum (-1)^{dim - ||cornerVector||_{1}} * I(x^{cornerVector})
    final DoubleType sum = new DoubleType();
    sum.setZero();
    // Convert from input to return type
    final Converter<I, DoubleType> conv = new RealDoubleConverter<>();
    final DoubleType valueAsDoubleType = new DoubleType();
    while (cursor.hasNext()) {
        final I value = cursor.next().copy();
        conv.convert(value, valueAsDoubleType);
        // Obtain the cursor position encoded as corner vector
        final int cornerInteger = cursor.getCornerRepresentation();
        // Determine if the value has to be added (factor==1) or subtracted
        // (factor==-1)
        final DoubleType factor = new DoubleType(Math.pow(-1.0d, dimensions - IntegralMean.norm(cornerInteger)));
        valueAsDoubleType.mul(factor);
        sum.add(valueAsDoubleType);
    }
    output.set(sum);
}

Example 64 with Cursor

use of net.imglib2.Cursor in project imagej-ops by imagej.

the class DefaultFrangi method frangi.

private final void frangi(RandomAccessibleInterval<T> in, RandomAccessibleInterval<U> out, int step) {
    // create denominators used for gaussians later.
    double ad = 2 * alpha * alpha;
    double bd = 2 * beta * beta;
    // OutOfBoundsMirrorStrategy for use when the cursor reaches the edges.
    OutOfBoundsMirrorFactory<T, RandomAccessibleInterval<T>> osmf = new OutOfBoundsMirrorFactory<T, RandomAccessibleInterval<T>>(Boundary.SINGLE);
    Cursor<T> cursor = Views.iterable(in).localizingCursor();
    Matrix hessian = new Matrix(in.numDimensions(), in.numDimensions());
    // use three RandomAccess<T> Objects to find the values needed to calculate
    // the
    // second derivatives.
    RandomAccess<T> current = osmf.create(in);
    RandomAccess<T> behind = osmf.create(in);
    RandomAccess<T> ahead = osmf.create(in);
    RandomAccess<U> outputRA = out.randomAccess();
    while (cursor.hasNext()) {
        cursor.fwd();
        // calculate the hessian
        for (int m = 0; m < in.numDimensions(); m++) {
            for (int n = 0; n < in.numDimensions(); n++) {
                current.setPosition(cursor);
                ahead.setPosition(cursor);
                behind.setPosition(cursor);
                // move one behind to take the first derivative
                behind.move(-step, m);
                if (m != n)
                    behind.move(-step, n);
                // take the derivative between the two points
                double derivativeA = derive(behind.get().getRealDouble(), current.get().getRealDouble(), getDistance(behind, current, in.numDimensions()));
                // move one ahead to take the other first derivative
                ahead.move(step, m);
                if (m != n)
                    ahead.move(step, n);
                // take the derivative between the two points
                double derivativeB = derive(current.get().getRealDouble(), ahead.get().getRealDouble(), getDistance(current, ahead, in.numDimensions()));
                // take the second derivative using the two first derivatives
                double derivative2 = derive(derivativeA, derivativeB, getDistance(behind, ahead, in.numDimensions()));
                hessian.set(m, n, derivative2);
            }
        }
        // find the FrobeniusNorm (used later)
        double s = hessian.normF();
        double cn = -(s * s);
        // find and sort the eigenvalues and eigenvectors of the Hessian
        EigenvalueDecomposition e = hessian.eig();
        double[] eigenvaluesArray = e.getRealEigenvalues();
        ArrayList<Double> eigenvaluesArrayList = new ArrayList<Double>();
        for (double d : eigenvaluesArray) eigenvaluesArrayList.add(d);
        Collections.sort(eigenvaluesArrayList, Comparator.comparingDouble(Math::abs));
        // vesselness value
        double v = 0;
        if (in.numDimensions() == 2) {
            double c = 15;
            double cd = 2 * c * c;
            // lambda values
            double l1 = eigenvaluesArrayList.get(0);
            double al1 = Math.abs(l1);
            double l2 = eigenvaluesArrayList.get(1);
            double al2 = Math.abs(l2);
            // Check to see if the point is on a tubular structure.
            if (l2 < 0) {
                // ratio Rb
                double rb = al1 / al2;
                // values for ease of final calculation
                double bn = -(rb * rb);
                v = Math.exp(bn / bd) * (1 - Math.exp(cn / cd));
            }
        } else if (in.numDimensions() == 3) {
            double c = 200;
            double cd = 2 * c * c;
            // lambda values
            double l1 = eigenvaluesArrayList.get(0);
            double al1 = Math.abs(l1);
            double l2 = eigenvaluesArrayList.get(1);
            double al2 = Math.abs(l2);
            double l3 = eigenvaluesArrayList.get(2);
            double al3 = Math.abs(l3);
            /*
				 * N.B. This conditional statement only takes into account the sign 
				 * on the third-smallest eigenvalue, not both the second and third as 
				 * described in the paper. Original versions of this filter took the
				 * signs of both into account, but the result of the filter was an 
				 * empty image. Only by removing the condition of the sign of the 
				 * second eigenvalue were we able to obtain results that matched 
				 * human expectations of a 3-D version of the filter. This conditional
				 * in particular achieved results best aligned with human expectation.
				 */
            if (l3 < 0) {
                // ratios Rb and Ra
                double rb = al1 / Math.sqrt(al2 * al3);
                double ra = al2 / al3;
                // values for ease of final calculation
                double an = -(ra * ra);
                double bn = -(rb * rb);
                // System.out.println(l1 + "  " + l2 + "  " + l3);
                v = (1 - Math.exp(an / ad)) * Math.exp(bn / bd) * (1 - Math.exp(cn / cd));
            }
        } else
            throw new RuntimeException("Currently only 2 or 3 dimensional images are supported");
        if (!Double.isNaN(v)) {
            maximumVesselness = Math.max(v, maximumVesselness);
            minimumVesselness = Math.min(v, minimumVesselness);
        }
        outputRA.setPosition(cursor);
        outputRA.get().setReal(v);
    }
}

Example 65 with Cursor

use of net.imglib2.Cursor in project imagej-ops by imagej.

the class DefaultContour method calculate.

@Override
public Polygon2D calculate(final RandomAccessibleInterval<B> input) {
    List<RealPoint> p = new ArrayList<>();
    final B var = Util.getTypeFromInterval(input).createVariable();
    final RandomAccess<B> raInput = Views.extendValue(input, var).randomAccess();
    final Cursor<B> cInput = Views.flatIterable(input).cursor();
    final ClockwiseMooreNeighborhoodIterator<B> cNeigh = new ClockwiseMooreNeighborhoodIterator<>(raInput);
    double[] position = new double[2];
    double[] startPos = new double[2];
    // find first true pixel
    while (cInput.hasNext()) {
        // we are looking for a true pixel
        if (cInput.next().get()) {
            raInput.setPosition(cInput);
            raInput.localize(startPos);
            // add to polygon
            p.add(new RealPoint(startPos[0], startPos[1]));
            // backtrack:
            raInput.move(-1, 0);
            cNeigh.reset();
            while (cNeigh.hasNext()) {
                if (cNeigh.next().get()) {
                    boolean specialBacktrack = false;
                    raInput.localize(position);
                    if (startPos[0] == position[0] && startPos[1] == position[1]) {
                        // startPoint was found.
                        if (useJacobs) {
                            // Jacobs stopping criteria
                            final int index = cNeigh.getIndex();
                            if (index == 1 || index == 0) {
                                // Jonathans refinement to
                                // non-terminating jacobs criteria
                                specialBacktrack = true;
                            } else if (index == 2 || index == 3) {
                                // way.
                                break;
                            }
                        // else criteria not fulfilled, continue.
                        } else {
                            break;
                        }
                    }
                    // add found point to polygon
                    p.add(new RealPoint(position[0], position[1]));
                    if (specialBacktrack) {
                        cNeigh.backtrackSpecial();
                    } else {
                        cNeigh.backtrack();
                    }
                }
            }
            // we only need to extract one contour.
            break;
        }
    }
    return new DefaultWritablePolygon2D(p);
}