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);
}
}
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);
}
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);
}
}
}
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);
}
}
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);
}
Aggregations