Search in sources :

Example 21 with Point

use of net.imglib2.Point 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);
    }
}
Also used : EigenvalueDecomposition(Jama.EigenvalueDecomposition) ArrayList(java.util.ArrayList) Matrix(Jama.Matrix) RandomAccessibleInterval(net.imglib2.RandomAccessibleInterval) OutOfBoundsMirrorFactory(net.imglib2.outofbounds.OutOfBoundsMirrorFactory)

Example 22 with Point

use of net.imglib2.Point 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);
}
Also used : RealPoint(net.imglib2.RealPoint) ArrayList(java.util.ArrayList) RealPoint(net.imglib2.RealPoint) DefaultWritablePolygon2D(net.imglib2.roi.geom.real.DefaultWritablePolygon2D)

Example 23 with Point

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

the class NonCirculantNormalizationFactor method createNormalizationImageSemiNonCirculant.

protected void createNormalizationImageSemiNonCirculant() {
    // k is the window size (valid image region)
    final int length = k.numDimensions();
    final long[] n = new long[length];
    final long[] nFFT = new long[length];
    // also referred to as object space size
    for (int d = 0; d < length; d++) {
        n[d] = k.dimension(d) + l.dimension(d) - 1;
    }
    for (int d = 0; d < length; d++) {
        nFFT[d] = imgConvolutionInterval.dimension(d);
    }
    FinalDimensions fd = new FinalDimensions(nFFT);
    // create the normalization image
    normalization = create.calculate(fd);
    // size of the measurement window
    final Point size = new Point(length);
    final long[] sizel = new long[length];
    for (int d = 0; d < length; d++) {
        size.setPosition(k.dimension(d), d);
        sizel[d] = k.dimension(d);
    }
    // starting point of the measurement window when it is centered in fft space
    final Point start = new Point(length);
    final long[] startl = new long[length];
    final long[] endl = new long[length];
    for (int d = 0; d < length; d++) {
        start.setPosition((nFFT[d] - k.dimension(d)) / 2, d);
        startl[d] = (nFFT[d] - k.dimension(d)) / 2;
        endl[d] = startl[d] + sizel[d] - 1;
    }
    // size of the object space
    final Point maskSize = new Point(length);
    final long[] maskSizel = new long[length];
    for (int d = 0; d < length; d++) {
        maskSize.setPosition(Math.min(n[d], nFFT[d]), d);
        maskSizel[d] = Math.min(n[d], nFFT[d]);
    }
    // starting point of the object space within the fft space
    final Point maskStart = new Point(length);
    final long[] maskStartl = new long[length];
    for (int d = 0; d < length; d++) {
        maskStart.setPosition((Math.max(0, nFFT[d] - n[d]) / 2), d);
        maskStartl[d] = (Math.max(0, nFFT[d] - n[d]) / 2);
    }
    final RandomAccessibleInterval<O> temp = Views.interval(normalization, new FinalInterval(startl, endl));
    final Cursor<O> normCursor = Views.iterable(temp).cursor();
    // draw a cube the size of the measurement space
    while (normCursor.hasNext()) {
        normCursor.fwd();
        normCursor.get().setReal(1.0);
    }
    final Img<O> tempImg = create.calculate(fd);
    // 3. correlate psf with the output of step 2.
    correlater.compute(normalization, tempImg);
    normalization = tempImg;
    final Cursor<O> cursorN = normalization.cursor();
    while (cursorN.hasNext()) {
        cursorN.fwd();
        if (cursorN.get().getRealFloat() <= 1e-3f) {
            cursorN.get().setReal(1.0f);
        }
    }
}
Also used : FinalDimensions(net.imglib2.FinalDimensions) FinalInterval(net.imglib2.FinalInterval) Point(net.imglib2.Point) Point(net.imglib2.Point)

Example 24 with Point

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

the class FFTTest method placeSphereInCenter.

/**
 * utility that places a sphere in the center of the image
 *
 * @param img
 */
private void placeSphereInCenter(final Img<FloatType> img) {
    final Point center = new Point(img.numDimensions());
    for (int d = 0; d < img.numDimensions(); d++) center.setPosition(img.dimension(d) / 2, d);
    final HyperSphere<FloatType> hyperSphere = new HyperSphere<>(img, center, 2);
    for (final FloatType value : hyperSphere) {
        value.setReal(1);
    }
}
Also used : HyperSphere(net.imglib2.algorithm.region.hypersphere.HyperSphere) Point(net.imglib2.Point) Point(net.imglib2.Point) FloatType(net.imglib2.type.numeric.real.FloatType) ComplexFloatType(net.imglib2.type.numeric.complex.ComplexFloatType)

Example 25 with Point

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

the class ConvolveTest method testCreateAndConvolvePoints.

/**
 * tests fft based convolve
 */
@Test
public void testCreateAndConvolvePoints() {
    final int xSize = 128;
    final int ySize = 128;
    final int zSize = 128;
    int[] size = new int[] { xSize, ySize, zSize };
    Img<DoubleType> phantom = ops.create().img(size);
    RandomAccess<DoubleType> randomAccess = phantom.randomAccess();
    randomAccess.setPosition(new long[] { xSize / 2, ySize / 2, zSize / 2 });
    randomAccess.get().setReal(255.0);
    randomAccess.setPosition(new long[] { xSize / 4, ySize / 4, zSize / 4 });
    randomAccess.get().setReal(255.0);
    Point location = new Point(phantom.numDimensions());
    location.setPosition(new long[] { 3 * xSize / 4, 3 * ySize / 4, 3 * zSize / 4 });
    HyperSphere<DoubleType> hyperSphere = new HyperSphere<>(phantom, location, 5);
    for (DoubleType value : hyperSphere) {
        value.setReal(16);
    }
    // create psf using the gaussian kernel op (alternatively PSF could be an
    // input to the script)
    RandomAccessibleInterval<DoubleType> psf = ops.create().kernelGauss(new double[] { 5, 5, 5 }, new DoubleType());
    // convolve psf with phantom
    RandomAccessibleInterval<DoubleType> convolved = ops.filter().convolve(phantom, psf);
    DoubleType sum = new DoubleType();
    DoubleType max = new DoubleType();
    DoubleType min = new DoubleType();
    ops.stats().sum(sum, Views.iterable(convolved));
    ops.stats().max(max, Views.iterable(convolved));
    ops.stats().min(min, Views.iterable(convolved));
    assertEquals(sum.getRealDouble(), 8750.00, 0.001);
    assertEquals(max.getRealDouble(), 3.155, 0.001);
    assertEquals(min.getRealDouble(), 2.978E-7, 0.001);
}
Also used : DoubleType(net.imglib2.type.numeric.real.DoubleType) HyperSphere(net.imglib2.algorithm.region.hypersphere.HyperSphere) Point(net.imglib2.Point) Point(net.imglib2.Point) AbstractOpTest(net.imagej.ops.AbstractOpTest) Test(org.junit.Test)

Aggregations

ArrayList (java.util.ArrayList)10 RealPoint (net.imglib2.RealPoint)10 Point (net.imglib2.Point)8 Test (org.junit.Test)7 RandomAccessibleInterval (net.imglib2.RandomAccessibleInterval)5 DoubleType (net.imglib2.type.numeric.real.DoubleType)5 Point (com.google.monitoring.v3.Point)4 AbstractFeatureTest (net.imagej.ops.features.AbstractFeatureTest)4 HyperSphere (net.imglib2.algorithm.region.hypersphere.HyperSphere)4 Polygon2D (net.imglib2.roi.geom.real.Polygon2D)4 FloatType (net.imglib2.type.numeric.real.FloatType)4 Metric (com.google.api.Metric)3 TimeInterval (com.google.monitoring.v3.TimeInterval)3 TimeSeries (com.google.monitoring.v3.TimeSeries)3 TypedValue (com.google.monitoring.v3.TypedValue)3 List (java.util.List)3 Interval (net.imglib2.Interval)3 MonitoredResource (com.google.api.MonitoredResource)2 MetricServiceClient (com.google.cloud.monitoring.v3.MetricServiceClient)2 CreateTimeSeriesRequest (com.google.monitoring.v3.CreateTimeSeriesRequest)2