use of net.imglib2.roi.labeling.LabelRegionCursor in project vcell by virtualcell.
the class ConstructTIRFGeometry method run.
@Override
public void run() {
// Calculate constant d in TIRF exponential decay function
// Angle of incidence in radians
theta = theta * 2 * Math.PI / 360;
// Refractive index of glass
final double n1 = 1.52;
// Refractive index of cytosol
final double n2 = 1.38;
final double d = lambda * Math.pow((Math.pow(n1, 2) * Math.pow(Math.sin(theta), 2) - Math.pow(n2, 2)), -0.5) / (4 * Math.PI);
System.out.println("d: " + d);
final double fluorPerMolecule = 250;
// Get frame of interest to define geometry
long maxX = data.dimension(0) - 1;
long maxY = data.dimension(1) - 1;
Interval interval = Intervals.createMinMax(0, 0, sliceIndex, maxX, maxY, sliceIndex);
RandomAccessibleInterval<T> croppedRAI = ops.transform().crop(data, interval, true);
// Subtract lowest pixel value
IterableInterval<T> dataII = Views.iterable(croppedRAI);
double min = ops.stats().min(dataII).getRealDouble();
Cursor<T> dataCursor = dataII.cursor();
while (dataCursor.hasNext()) {
double val = dataCursor.next().getRealDouble();
dataCursor.get().setReal(val - min);
}
// Perform Gaussian blur
RandomAccessibleInterval<T> blurredRAI = ops.filter().gauss(croppedRAI, 2);
IterableInterval<T> blurredII = Views.iterable(blurredRAI);
// Segment slice by threshold and fill holes
IterableInterval<BitType> thresholded = ops.threshold().huang(blurredII);
Img<BitType> thresholdedImg = ops.convert().bit(thresholded);
RandomAccessibleInterval<BitType> thresholdedRAI = ops.morphology().fillHoles(thresholdedImg);
// Get the largest region
RandomAccessibleInterval<LabelingType<ByteType>> labeling = ops.labeling().cca(thresholdedRAI, ConnectedComponents.StructuringElement.EIGHT_CONNECTED);
LabelRegions<ByteType> labelRegions = new LabelRegions<>(labeling);
Iterator<LabelRegion<ByteType>> iterator = labelRegions.iterator();
LabelRegion<ByteType> maxRegion = iterator.next();
while (iterator.hasNext()) {
LabelRegion<ByteType> currRegion = iterator.next();
if (currRegion.size() > maxRegion.size()) {
maxRegion = currRegion;
}
}
// Generate z index map
double iMax = ops.stats().max(dataII).getRealDouble();
Img<UnsignedShortType> dataImg = ops.convert().uint16(dataII);
Img<UnsignedShortType> zMap = ops.convert().uint16(ops.create().img(dataII));
LabelRegionCursor cursor = maxRegion.localizingCursor();
RandomAccess<UnsignedShortType> zMapRA = zMap.randomAccess();
RandomAccess<UnsignedShortType> dataRA = dataImg.randomAccess();
while (cursor.hasNext()) {
cursor.fwd();
zMapRA.setPosition(cursor);
dataRA.setPosition(cursor);
double val = dataRA.get().getRealDouble();
// Log of 0 is undefined
if (val < 1) {
val = 1;
}
int z = (int) Math.round(-d * Math.log(val / iMax) / zRes);
zMapRA.get().set(z);
}
System.out.println("6");
// Use map to construct 3D geometry
// Add 5 slices of padding on top
int maxZ = (int) ops.stats().max(zMap).getRealDouble() + 5;
long[] resultDimensions = { maxX + 1, maxY + 1, maxZ };
Img<BitType> result = new ArrayImgFactory<BitType>().create(resultDimensions, new BitType());
RandomAccess<BitType> resultRA = result.randomAccess();
System.out.println(maxZ);
cursor.reset();
while (cursor.hasNext()) {
cursor.fwd();
zMapRA.setPosition(cursor);
int zIndex = zMapRA.get().get();
int[] position = { cursor.getIntPosition(0), cursor.getIntPosition(1), zIndex };
while (position[2] < maxZ) {
resultRA.setPosition(position);
resultRA.get().set(true);
position[2]++;
}
}
output = datasetService.create(result);
CalibratedAxis[] axes = new DefaultLinearAxis[] { new DefaultLinearAxis(Axes.X), new DefaultLinearAxis(Axes.Y), new DefaultLinearAxis(Axes.Z) };
output.setAxes(axes);
System.out.println("Done constructing geometry");
}
Aggregations