use of net.imglib2.FinalInterval in project imagej-ops by imagej.
the class DefaultVoxelization3D method calculate.
@Override
public RandomAccessibleInterval<BitType> calculate(Mesh input) {
Img<BitType> outImg = ops.create().img(new FinalInterval(width, height, depth), new BitType());
DefaultMesh dMesh = (DefaultMesh) input;
Set<RealLocalizable> verts = dMesh.getVertices();
RealPoint minPoint = new RealPoint(verts.iterator().next());
RealPoint maxPoint = new RealPoint(verts.iterator().next());
for (RealLocalizable v : verts) {
if (v.getDoublePosition(0) < minPoint.getDoublePosition(0))
minPoint.setPosition(v.getDoublePosition(0), 0);
if (v.getDoublePosition(1) < minPoint.getDoublePosition(1))
minPoint.setPosition(v.getDoublePosition(1), 1);
if (v.getDoublePosition(2) < minPoint.getDoublePosition(2))
minPoint.setPosition(v.getDoublePosition(2), 2);
if (v.getDoublePosition(0) > maxPoint.getDoublePosition(0))
maxPoint.setPosition(v.getDoublePosition(0), 0);
if (v.getDoublePosition(1) > maxPoint.getDoublePosition(1))
maxPoint.setPosition(v.getDoublePosition(1), 1);
if (v.getDoublePosition(2) > maxPoint.getDoublePosition(2))
maxPoint.setPosition(v.getDoublePosition(2), 2);
}
RealPoint dimPoint = new RealPoint((maxPoint.getDoublePosition(0) - minPoint.getDoublePosition(0)), (maxPoint.getDoublePosition(1) - minPoint.getDoublePosition(1)), (maxPoint.getDoublePosition(2) - minPoint.getDoublePosition(2)));
double[] stepSizes = new double[3];
stepSizes[0] = dimPoint.getDoublePosition(0) / width;
stepSizes[1] = dimPoint.getDoublePosition(1) / height;
stepSizes[2] = dimPoint.getDoublePosition(2) / depth;
double[] voxelHalfsize = new double[3];
for (int k = 0; k < stepSizes.length; k++) voxelHalfsize[k] = stepSizes[k] / 2.0;
for (Facet f : dMesh.getFacets()) {
TriangularFacet tri = (TriangularFacet) f;
Vector3D v1 = tri.getP0();
Vector3D v2 = tri.getP1();
Vector3D v3 = tri.getP2();
double[] minSubBoundary = new double[] { Math.min(Math.min(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), Math.min(Math.min(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), Math.min(Math.min(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) };
double[] maxSubBoundary = new double[] { Math.max(Math.max(v1.getX(), v2.getX()), v3.getX()) - minPoint.getDoublePosition(0), Math.max(Math.max(v1.getY(), v2.getY()), v3.getY()) - minPoint.getDoublePosition(1), Math.max(Math.max(v1.getZ(), v2.getZ()), v3.getZ()) - minPoint.getDoublePosition(2) };
// Should use the
RandomAccess<BitType> ra = outImg.randomAccess();
// interval
// implementation
// for speed
long[] indices = new long[3];
for (indices[0] = (long) Math.floor(minSubBoundary[0] / stepSizes[0]); indices[0] < Math.floor(maxSubBoundary[0] / stepSizes[0]); indices[0]++) {
for (indices[1] = (long) Math.floor(minSubBoundary[1] / stepSizes[1]); indices[1] < Math.floor(maxSubBoundary[1] / stepSizes[1]); indices[1]++) {
for (indices[2] = (long) Math.floor(minSubBoundary[2] / stepSizes[2]); indices[2] < Math.floor(maxSubBoundary[2] / stepSizes[2]); indices[2]++) {
ra.setPosition(indices);
if (// Don't check if voxel is already
!ra.get().get()) // filled
{
double[] voxelCenter = new double[3];
for (int k = 0; k < 3; k++) voxelCenter[k] = indices[k] * stepSizes[k] + voxelHalfsize[k];
if (triBoxOverlap(voxelCenter, voxelHalfsize, v1, v2, v3) == 1) {
ra.get().set(true);
}
}
}
}
}
}
return outImg;
}
use of net.imglib2.FinalInterval in project imagej-ops by imagej.
the class DefaultCreateKernelGabor method calculate.
@Override
public RandomAccessibleInterval<T> calculate(final double[] sigmas, final double[] period) {
// both input arrays must be of the same length
if (sigmas.length != period.length)
throw new IllegalArgumentException("Params length mismatch: The number " + "of sigmas must match the dimensionality of the period vector.");
// NB: sigma==0 indicates no filtering along its axis
for (final double s : sigmas) if (s < 0.0)
throw new IllegalArgumentException("Input sigma must be non-negative.");
// the size and center of the output image
final long[] dims = new long[sigmas.length];
final long[] centre = new long[sigmas.length];
for (int d = 0; d < dims.length; d++) {
dims[d] = Math.max(3, (2 * (int) (3 * sigmas[d] + 0.5) + 1));
centre[d] = (int) (dims[d] / 2);
}
// prepare the output image
final RandomAccessibleInterval<T> out = createImgOp.calculate(new FinalInterval(dims));
// calculate the squared length of the period vector
double perLengthSq = 0.0;
for (int d = 0; d < period.length; d++) perLengthSq += period[d] * period[d];
// fill the output image
final Cursor<T> cursor = Views.iterable(out).cursor();
while (cursor.hasNext()) {
cursor.fwd();
// obtain the current coordinate (use dims to store it)
cursor.localize(dims);
// to calculate current Gabor kernel value
double GaussExp = 0.0;
double freqPart = 0.0;
// but produce no Gaussian envelope for axes for which sigma==0
// no blocking by default
double blockingExp = 1.0;
// sweep over all dimensions to determine voxel value
for (int d = 0; d < dims.length; d++) {
final double dx = dims[d] - centre[d];
if (sigmas[d] > 0.)
// normal case: cummulate exp's argument
GaussExp += (dx * dx) / (sigmas[d] * sigmas[d]);
else if (dx != 0.)
// sigmas[d] == 0 && we are off the blocking axis
blockingExp = 0.f;
// cummulates scalar product...
freqPart += dx * period[d];
}
GaussExp = Math.exp(-0.5 * GaussExp) * blockingExp;
freqPart = 6.28318 * freqPart / perLengthSq;
// compose the real value finally
cursor.get().setReal(GaussExp * Math.cos(freqPart));
// TODO NB: is it faster to determine type or calculate the math (possible uselessly)
if (!(typeVar instanceof RealType<?>))
// set then the imaginary part of the kernel too
cursor.get().setImaginary(GaussExp * Math.sin(freqPart));
}
return out;
}
use of net.imglib2.FinalInterval in project imagej-ops by imagej.
the class DefaultCreateKernelGauss method calculate.
@Override
public RandomAccessibleInterval<T> calculate(double[] input) {
final double[] sigmaPixels = new double[input.length];
final long[] dims = new long[input.length];
final double[][] kernelArrays = new double[input.length][];
for (int d = 0; d < input.length; d++) {
sigmaPixels[d] = input[d];
dims[d] = Math.max(3, (2 * (int) (3 * sigmaPixels[d] + 0.5) + 1));
kernelArrays[d] = Util.createGaussianKernel1DDouble(sigmaPixels[d], true);
}
final RandomAccessibleInterval<T> out = createOp.calculate(new FinalInterval(dims));
final Cursor<T> cursor = Views.iterable(out).cursor();
while (cursor.hasNext()) {
cursor.fwd();
double result = 1.0f;
for (int d = 0; d < input.length; d++) {
result *= kernelArrays[d][cursor.getIntPosition(d)];
}
cursor.get().setReal(result);
}
return out;
}
use of net.imglib2.FinalInterval in project imagej-ops by imagej.
the class CreateKernelSobel method calculate.
@Override
public RandomAccessibleInterval<T> calculate() {
long[] dim = new long[4];
dim[0] = 3;
dim[1] = 1;
for (int k = 2; k < dim.length; k++) {
dim[k] = 1;
}
dim[dim.length - 1] = 2;
RandomAccessibleInterval<T> output = createOp.calculate(new FinalInterval(dim));
final Cursor<T> cursor = Views.iterable(output).cursor();
int i = 0;
while (cursor.hasNext()) {
cursor.fwd();
cursor.get().setReal(values[i]);
i++;
}
return output;
}
use of net.imglib2.FinalInterval in project imagej-ops by imagej.
the class DefaultCreateKernelBiGauss method calculate.
@Override
public RandomAccessibleInterval<T> calculate(final double[] sigmas, final Integer dimensionality) {
// both sigmas must be available
if (sigmas.length < 2)
throw new IllegalArgumentException("Two sigmas (for inner and outer Gauss)" + " must be supplied.");
// both sigmas must be reasonable
if (sigmas[0] <= 0 || sigmas[1] <= 0)
throw new IllegalArgumentException("Input sigmas must be both positive.");
// dimension as well...
if (dimensionality <= 0)
throw new IllegalArgumentException("Input dimensionality must both positive.");
// the size and center of the output image
final long[] dims = new long[dimensionality];
final long[] centre = new long[dimensionality];
// time-saver... (must hold now: dimensionality > 0)
dims[0] = Math.max(3, (2 * (int) (sigmas[0] + 2 * sigmas[1] + 0.5) + 1));
centre[0] = (int) (dims[0] / 2);
// fill the size and center arrays
for (int d = 1; d < dims.length; d++) {
dims[d] = dims[0];
centre[d] = centre[0];
}
// prepare some scaling constants
// eq. (6)
final double k = (sigmas[1] / sigmas[0]) * (sigmas[1] / sigmas[0]);
// eq. (9)
final double c0 = 0.24197 * ((sigmas[1] / sigmas[0]) - 1.0) / sigmas[0];
// 0.24197 = 1/sqrt(2*PI*e) = 1/sqrt(2*PI) * exp(-0.5)
final double[] C = { 1.0 / (2.50663 * sigmas[0]), 1.0 / (2.50663 * sigmas[1]) };
// 2.50663 = sqrt(2*PI)
// prepare squared input sigmas
final double[] sigmasSq = { sigmas[0] * sigmas[0], sigmas[1] * sigmas[1] };
// prepare the output image
final RandomAccessibleInterval<T> out = createImgOp.calculate(new FinalInterval(dims));
// fill the output image
final Cursor<T> cursor = Views.iterable(out).cursor();
while (cursor.hasNext()) {
cursor.fwd();
// obtain the current coordinate (use dims to store it)
cursor.localize(dims);
// calculate distance from the image centre
// TODO: can JVM reuse this var or is it allocated again and again (and multipling in the memory)?
double dist = 0.;
for (int d = 0; d < dims.length; d++) {
final double dx = dims[d] - centre[d];
dist += dx * dx;
}
// dist = Math.sqrt(dist); -- gonna work with squared distance
// which of the two Gaussians should we use?
double val = 0.;
if (dist < sigmasSq[0]) {
// the inner one
val = C[0] * Math.exp(-0.5 * dist / sigmasSq[0]) + c0;
} else {
// the outer one, get new distance first:
dist = Math.sqrt(dist) - (sigmas[0] - sigmas[1]);
dist *= dist;
val = k * C[1] * Math.exp(-0.5 * dist / sigmasSq[1]);
}
// compose the real value finally
cursor.get().setReal(val);
}
return out;
}
Aggregations