Examples with OutOfBoundsMirrorFactory net.imglib2.outofbounds.OutOfBoundsMirrorFactory used on opensource projects

Example 36 with OutOfBoundsMirrorFactory

use of net.imglib2.outofbounds.OutOfBoundsMirrorFactory 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);
    }
}