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

Example 16 with Point

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

the class DefaultDetectJunctionsTest method TestRegression.

@Test
public void TestRegression() {
    List<DefaultWritablePolyline> lines = new ArrayList<>();
    double threshold = Math.sqrt(2);
    // create first polyline (horizontal)
    List<RealPoint> list1 = new ArrayList<>();
    RealPoint p = new RealPoint(-5, 0);
    for (int i = 0; i < 10; i++) {
        p.move(1, 0);
        list1.add(new RealPoint(p));
    }
    lines.add(new DefaultWritablePolyline(list1));
    // create second polyline (vertical)
    List<RealPoint> list2 = new ArrayList<>();
    p.setPosition(0, 0);
    p.setPosition(-5, 1);
    for (int i = 0; i < 10; i++) {
        p.move(1, 1);
        list2.add(new RealPoint(p));
    }
    lines.add(new DefaultWritablePolyline(list2));
    // create third polyline (diagonal)
    List<RealPoint> list3 = new ArrayList<>();
    p.setPosition(-15, 0);
    p.setPosition(-15, 1);
    for (int i = 0; i < 20; i++) {
        p.move(1, 0);
        p.move(1, 1);
        list3.add(new RealPoint(p));
    }
    lines.add(new DefaultWritablePolyline(list3));
    // create fourth polyline (vertical, different intersection point)
    List<RealPoint> list4 = new ArrayList<>();
    p.setPosition(-11, 0);
    p.setPosition(-18, 1);
    for (int i = 0; i < 7; i++) {
        p.move(1, 1);
        list4.add(new RealPoint(p));
    }
    lines.add(new DefaultWritablePolyline(list4));
    List<RealPoint> results;
    results = (List<RealPoint>) ops.run(net.imagej.ops.segment.detectJunctions.DefaultDetectJunctions.class, lines, threshold);
    List<RealPoint> expected = new ArrayList<>();
    expected.add(new RealPoint(0, 0));
    expected.add(new RealPoint(-11, -11));
    for (int i = 0; i < results.size(); i++) {
        assertEquals(results.get(i).getDoublePosition(0), expected.get(i).getDoublePosition(0), 0);
        assertEquals(results.get(i).getDoublePosition(1), expected.get(i).getDoublePosition(1), 0);
    }
}

Example 17 with Point

use of net.imglib2.Point 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 18 with Point

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

the class DefaultDetectJunctions method calculate.

@Override
public List<RealPoint> calculate(List<? extends WritablePolyline> input) {
    // output that allows for both split polyline inputs and a
    // realPointCollection for our junctions.
    List<RealPoint> output = new ArrayList<>();
    for (int first = 0; first < input.size() - 1; first++) {
        WritablePolyline firstLine = input.get(first);
        for (int second = first + 1; second < input.size(); second++) {
            WritablePolyline secondLine = input.get(second);
            // interval containing both plines
            Interval intersect = Intervals.intersect(slightlyEnlarge(firstLine, 2), slightlyEnlarge(secondLine, 2));
            // each other.
            if (Intervals.isEmpty(intersect))
                continue;
            // create an arraylist to contain all of the junctions for these two
            // lines (so that we can filter the junctions before putting them in
            // output).
            ArrayList<RealPoint> currentPairJunctions = new ArrayList<>();
            for (int p = 0; p < firstLine.numVertices() - 1; p++) {
                for (int q = 0; q < secondLine.numVertices() - 1; q++) {
                    RealLocalizableRealPositionable p1 = firstLine.vertex(p);
                    RealLocalizableRealPositionable p2 = firstLine.vertex(p + 1);
                    RealLocalizableRealPositionable q1 = secondLine.vertex(q);
                    RealLocalizableRealPositionable q2 = secondLine.vertex(q + 1);
                    // special cases if both lines are vertical
                    boolean pVertical = Math.round(p1.getDoublePosition(0)) == Math.round(p2.getDoublePosition(0));
                    boolean qVertical = Math.round(q1.getDoublePosition(0)) == Math.round(q2.getDoublePosition(0));
                    // intersection point between the lines created by line segments p
                    // and q.
                    double[] intersectionPoint = new double[2];
                    // since they are parallel and cannot be the same.
                    if (pVertical && qVertical) {
                        parallelRoutine(p1, p2, q1, q2, currentPairJunctions, true);
                        continue;
                    } else if (pVertical) {
                        double mq = (q2.getDoublePosition(1) - q1.getDoublePosition(1)) / (q2.getDoublePosition(0) - q1.getDoublePosition(0));
                        double bq = (q1.getDoublePosition(1) - mq * q1.getDoublePosition(0));
                        double x = p1.getDoublePosition(0);
                        double y = mq * x + bq;
                        intersectionPoint[0] = x;
                        intersectionPoint[1] = y;
                    } else if (qVertical) {
                        double mp = (p2.getDoublePosition(1) - p1.getDoublePosition(1)) / (p2.getDoublePosition(0) - p1.getDoublePosition(0));
                        double bp = (p1.getDoublePosition(1) - mp * p1.getDoublePosition(0));
                        double x = q1.getDoublePosition(0);
                        double y = mp * x + bp;
                        intersectionPoint[0] = x;
                        intersectionPoint[1] = y;
                    } else {
                        double mp = (p2.getDoublePosition(1) - p1.getDoublePosition(1)) / (p2.getDoublePosition(0) - p1.getDoublePosition(0));
                        double mq = (q2.getDoublePosition(1) - q1.getDoublePosition(1)) / (q2.getDoublePosition(0) - q1.getDoublePosition(0));
                        if (mp == mq) {
                            parallelRoutine(p1, p2, q1, q2, currentPairJunctions, false);
                            continue;
                        }
                        double bp = (p2.getDoublePosition(1) - mp * p2.getDoublePosition(0));
                        double bq = (q2.getDoublePosition(1) - mq * q2.getDoublePosition(0));
                        // point of intersection of lines created by line segments p and
                        // q.
                        double x = (bq - bp) / (mp - mq);
                        double y = mp * x + bp;
                        intersectionPoint[0] = x;
                        intersectionPoint[1] = y;
                    }
                    // find the distance from the intersection point to both line
                    // segments, and the length of the line segments.
                    double distp1 = getDistance(intersectionPoint, p1);
                    double distp2 = getDistance(intersectionPoint, p2);
                    double distq1 = getDistance(intersectionPoint, q1);
                    double distq2 = getDistance(intersectionPoint, q2);
                    // max distance from line segment to intersection point
                    double maxDist = Math.max(Math.min(distp1, distp2), Math.min(distq1, distq2));
                    // these lines are close enough to form a junction
                    if (maxDist <= threshold)
                        currentPairJunctions.add(new RealPoint(intersectionPoint));
                }
            }
            // filter out the current pair's junctions by removing duplicates and
            // then averaging all remaining nearby junctions
            filterJunctions(currentPairJunctions);
            // add the filtered junctions to the output list.
            for (RealPoint point : currentPairJunctions) output.add(point);
        }
    }
    // filter the junctions -- for each set of junctions that seem vaguely
    // similar, pick out the best one-
    filterJunctions(output);
    return output;
}

Example 19 with Point

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

the class DefaultDetectRidges method getNextPoint.

/**
 * Recursively determines the next line point and adds it to the running list
 * of line points.
 *
 * @param gradientRA - the {@link RandomAccess} of the gradient image.
 * @param pRA - the {@link RandomAccess} of the eigenvector image.
 * @param nRA - the {@link RandomAccess} of the subpixel line location image.
 * @param points - the {@link ArrayList} containing the line points.
 * @param octant - integer denoting the octant of the last gradient vector,
 *          oriented with 1 being 0 degrees and increasing in the
 *          counterclockwise direction.
 * @param lastnx - the x component of the gradient vector of the last line
 *          point.
 * @param lastny - the y component of the gradient vector of the last line
 *          point.
 * @param lastpx - the x component of the subpixel line location of the last
 *          line point.
 * @param lastpy - the y component of the subpixel line location of the last
 *          line point.
 */
private void getNextPoint(RandomAccess<DoubleType> gradientRA, RandomAccess<DoubleType> pRA, RandomAccess<DoubleType> nRA, List<RealPoint> points, int octant, double lastnx, double lastny, double lastpx, double lastpy) {
    Point currentPos = new Point(gradientRA);
    // variables for the best line point of the three.
    Point salientPoint = new Point(gradientRA);
    double salientnx = 0;
    double salientny = 0;
    double salientpx = 0;
    double salientpy = 0;
    double bestSalience = Double.MAX_VALUE;
    boolean lastPointInLine = true;
    // check the three possible points that could continue the line, starting at
    // the octant after the given octant and rotating clockwise around the
    // current pixel.
    double lastAngle = RidgeDetectionUtils.getAngle(lastnx, lastny);
    for (int i = 1; i < 4; i++) {
        int[] modifier = RidgeDetectionUtils.getOctantCoords(octant + i);
        gradientRA.move(modifier[0], 0);
        gradientRA.move(modifier[1], 1);
        // there.
        if (gradientRA.get().get() > lowerThreshold) /*&& isMaxRA.get().get() > 0*/
        {
            long[] vectorArr = { gradientRA.getLongPosition(0), gradientRA.getLongPosition(1), 0 };
            nRA.setPosition(vectorArr);
            double nx = nRA.get().get();
            nRA.fwd(2);
            double ny = nRA.get().get();
            pRA.setPosition(vectorArr);
            double px = pRA.get().get();
            pRA.fwd(2);
            double py = pRA.get().get();
            double currentAngle = RidgeDetectionUtils.getAngle(nx, ny);
            double subpixelDiff = Math.sqrt(Math.pow(px - lastpx, 2) + Math.pow(py - lastpy, 2));
            double angleDiff = Math.abs(currentAngle - lastAngle);
            lastPointInLine = false;
            // numbers relative to other potential line points.
            if (subpixelDiff + angleDiff < bestSalience) {
                // record the values of the new most salient pixel
                salientPoint = new Point(gradientRA);
                salientnx = nx;
                salientny = ny;
                salientpx = px;
                salientpy = py;
                bestSalience = subpixelDiff + angleDiff;
            }
            // set the values to zero so that they are not added to another line.
            gradientRA.get().set(0);
        }
        // reset our randomAccess for the next check
        gradientRA.setPosition(currentPos);
    }
    // set the current pixel to 0 in the first slice of eigenRA!
    gradientRA.get().setReal(0);
    // find the next line point as long as there is one to find
    if (!lastPointInLine) {
        // take the most salient point
        gradientRA.setPosition(salientPoint);
        points.add(RidgeDetectionUtils.get2DRealPoint(gradientRA.getDoublePosition(0) + salientpx, gradientRA.getDoublePosition(1) + salientpy));
        // the gradient vector itself refers to the greatest change in intensity,
        // and for a pixel on a line this vector will be perpendicular to the
        // direction of the line. But this vector can point to either the left or
        // the right of the line from the perspective of the detector, and there
        // is no guarantee that the vectors at line point will point off the same
        // side of the line. So if they point off different sides, set the current
        // vector by 180 degrees for the purposes of this detector. We set the
        // threshold for angle fixing just above 90 degrees since any lower would
        // prevent ridges curving.
        double potentialGradient = RidgeDetectionUtils.getAngle(salientnx, salientny);
        // even though they are close enough to satisfy.
        if (lastAngle < angleThreshold)
            lastAngle += 360;
        if (potentialGradient < angleThreshold)
            potentialGradient += 360;
        if (Math.abs(potentialGradient - lastAngle) > angleThreshold) {
            salientnx = -salientnx;
            salientny = -salientny;
        }
        // perform the operation again on the new end of the line being formed.
        getNextPoint(gradientRA, pRA, nRA, points, RidgeDetectionUtils.getOctant(salientnx, salientny), salientnx, salientny, salientpx, salientpy);
    }
}

Example 20 with Point

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

the class DefaultDetectRidges method calculate.

@Override
public List<? extends WritablePolyline> calculate(RandomAccessibleInterval<T> input) {
    double sigma = (width / (2 * Math.sqrt(3)));
    // generate the metadata images
    RidgeDetectionMetadata ridgeDetectionMetadata = new RidgeDetectionMetadata(input, sigma, lowerThreshold, higherThreshold);
    // retrieve the metadata images
    Img<DoubleType> p_values = ridgeDetectionMetadata.getPValues();
    Img<DoubleType> n_values = ridgeDetectionMetadata.getNValues();
    Img<DoubleType> gradients = ridgeDetectionMetadata.getGradients();
    // create RandomAccesses for the metadata images
    OutOfBoundsConstantValueFactory<DoubleType, RandomAccessibleInterval<DoubleType>> oscvf = new OutOfBoundsConstantValueFactory<>(new DoubleType(0));
    RandomAccess<DoubleType> pRA = oscvf.create(p_values);
    RandomAccess<DoubleType> nRA = oscvf.create(n_values);
    RandomAccess<DoubleType> gradientRA = oscvf.create(gradients);
    // create the output polyline list.
    List<DefaultWritablePolyline> lines = new ArrayList<>();
    // start at the point of greatest maximum absolute value
    gradientRA.setPosition(RidgeDetectionUtils.getMaxCoords(gradients, true));
    // loop through the maximum values of the image
    while (Math.abs(gradientRA.get().get()) > higherThreshold) {
        // create the List of points that will be used to make the polyline
        List<RealPoint> points = new ArrayList<>();
        // get all of the necessary metadata from the image.
        long[] eigenvectorPos = { gradientRA.getLongPosition(0), gradientRA.getLongPosition(1), 0 };
        // obtain the n-values
        nRA.setPosition(eigenvectorPos);
        double eigenx = nRA.get().getRealDouble();
        nRA.fwd(2);
        double eigeny = nRA.get().getRealDouble();
        // obtain the p-values
        pRA.setPosition(eigenvectorPos);
        double px = pRA.get().getRealDouble();
        pRA.fwd(2);
        double py = pRA.get().getRealDouble();
        // start the list by adding the current point, which is the most line-like
        // point on the polyline
        points.add(RidgeDetectionUtils.get2DRealPoint(gradientRA.getDoublePosition(0) + px, gradientRA.getDoublePosition(1) + py));
        // go in the direction to the left of the perpendicular value
        getNextPoint(gradientRA, pRA, nRA, points, RidgeDetectionUtils.getOctant(eigenx, eigeny), eigenx, eigeny, px, py);
        // flip the array list around so that we get one cohesive line
        gradientRA.setPosition(new long[] { eigenvectorPos[0], eigenvectorPos[1] });
        Collections.reverse(points);
        // go in the opposite direction as before.
        eigenx = -eigenx;
        eigeny = -eigeny;
        getNextPoint(gradientRA, pRA, nRA, points, RidgeDetectionUtils.getOctant(eigenx, eigeny), eigenx, eigeny, px, py);
        // set the value to 0 so that it is not reused.
        gradientRA.get().setReal(0);
        // list has fewer vertices than the parameter, then we do not report it.
        if (points.size() > ridgeLengthMin) {
            DefaultWritablePolyline pline = new DefaultWritablePolyline(points);
            lines.add(pline);
        }
        // find the next max absolute value
        gradientRA.setPosition(RidgeDetectionUtils.getMaxCoords(gradients, true));
    }
    return lines;
}