public void traceWalls() {
int width = control.getStack().getWidth();
int height = control.getStack().getHeight();
int depth = control.getStack().getDepth();
control.clearChannels();
SliceWriter lumenlog = new SliceWriter("lumens.raw", width, height, depth, 5);
control.addChannel(lumenlog, "Lumens");
SliceWriter walls = new SliceWriter("walls.raw", width, height, depth, 5);
control.addChannel(walls, "Walls");
SliceWriter skeletonLog = new SliceWriter("skeleton.raw", width, height, depth, 5);
control.addChannel(skeletonLog, "Skeleton");
keepTracking = true;
for (int z = 0; z < depth && keepTracking; z++) {
long time = System.currentTimeMillis();
short[] pixels = control.getSlice(z).clone();
short[] map = ImageFunctions.distance2d(pixels, width, height);
short[] mark = new short[width * height];
short[] sout = new short[pixels.length];
ArrayList<Lumen> lumens = new ArrayList<Lumen>();
for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
int index = x + y * width;
if (mark[index] == 0 && pixels[index] != 0 && !isReserved(pixels[index])) {
short[] newLumenMask = new short[pixels.length];
short fid = pixels[index];
int count = 0;
for (int i = 0; i < newLumenMask.length; i++) {
if (pixels[i] == fid) {
newLumenMask[i] = fid;
mark[i] = (short) 1;
count++;
}
}
Rectangle bounds = ImageFunctions.getBounds(newLumenMask, width, height, 12);
int regionID = 1;
short[] mark2 = new short[width * height];
ArrayList<FiberSlice> liste = new ArrayList<FiberSlice>();
for (int xx = 0; xx < width; xx++) {
for (int yy = 0; yy < height; yy++) {
if (mark2[xx + yy * width] == 0 && newLumenMask[xx + yy * width] == fid) {
regionID++;
ImageFunctions.floodFill2D(xx, yy, width, height, newLumenMask, mark2, (short) regionID);
short[] lumenPart = new short[newLumenMask.length];
int size = ImageFunctions.floodFill2D(xx, yy, width, height, newLumenMask, lumenPart, (short) regionID);
FiberSlice storedLumenPart = new FiberSlice(lumenPart, size, z, fid);
liste.add(storedLumenPart);
}
}
}
// Find path between the two largest lumen parts
if (liste.size() > 1) {
Collections.sort(liste);
// Get the two largest lumen parts
FiberSlice lumen1 = liste.remove(0);
FiberSlice lumen2 = liste.remove(0);
// find border of lumen 2
short[] lumenBorder = ImageFunctions.dilate(lumen2.lumendata.clone(), width, height);
lumenBorder = ImageFunctions.mask(lumenBorder, ImageFunctions.invertMask(lumen2.lumendata, true), true);
// find geodesic distance from lumen 1
short[] geodist = ImageFunctions.geodist(lumen1.lumendata, pixels, width, height);
for (int i = 0; i < geodist.length; i++) {
if (geodist[i] == 0) {
geodist[i] = 255;
}
}
// find path that connects lumen1 with lumen2 and
// goes through non-void part in the image
DistancePath bestPath = new DistancePath();
short[] path = bestPath.getPath(geodist,
lumenBorder, map, width, height);
for (int i = 0; i < path.length; i++) {
if (path[i] != 0) {
mark[i] = fid; // MUST AVOID FINDING NEW PIXELS AS NEW REGION
newLumenMask[i] = fid;
pixels[i] = fid;
}
}
}
// Calculate and store distance from lumen
short[] dmap = ImageFunctions.distance2d(ImageFunctions.invertMask(newLumenMask, false), width, height, bounds);
Lumen newLumen = new Lumen(newLumenMask, dmap, bounds, count, fid);
lumens.add(newLumen);
}
}
}
// Copy lumens into a separate file (only used for visual inspection
// later)
short[] lumenOutput = new short[pixels.length];
for (int i = 0; i < lumenOutput.length; i++) {
// copy only non-white pixels to output
if (!isReserved((short) i)) {
lumenOutput[i] = pixels[i];
}
}
lumenlog.putSlice(lumenOutput, z);
short[] waterShed = getWatershed(z);
for (Lumen lumen : lumens) {
Rectangle bounds = lumen.getBounds();
int startX = bounds.x;
int stopX = bounds.width + bounds.x;
int startY = bounds.y;
int stopY = bounds.height + bounds.y;
short[] currentWall = waterShed.clone();
short[] searchMask = ImageFunctions.dilate(lumen.getLumenMask(), width, height);
for (int q = 0; q < pixels.length; q++) {
if (lumen.getLumenMask()[q] != 0) {
searchMask[q] = 0;
}
}
// Join regions inside lumens
for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
int index = x + y * width;
if (currentWall[index] == 0 && searchMask[index] != 0) {
ImageFunctions.floodFill2D(currentWall.clone(), currentWall, x, y, width, height, (short) 1);
}
}
}
// Fill in watershed borders
for (int x = 0; x < pixels.length; x++) {
if (currentWall[x] == (short) 255) {
currentWall[x] = 0;
}
}
// dilate to fill borders
for (int i = 0; i < 3; i++) {
currentWall = ImageFunctions.dilate(currentWall.clone(), currentWall, width, height, lumen.getBounds());
}
// Obtain the original pixels ( without connecting paths written to them)
pixels = control.getSlice(z);
// get median wallsize
short ff = (short) 0xff;
ArrayList<Integer> sizeList = new ArrayList<Integer>();
for (int x = startX; x < stopX; x++) {
for (int y = startY; y < stopY; y++) {
int i = x + y * width;
if (pixels[i] == ff && lumen.getDistanceMap()[i] < 9) {
sizeList.add((int) map[i]);
}
}
}
Collections.sort(sizeList);
if (sizeList.size() == 0) {
System.err.println("warning sizelist empty!");
break;
}
int medianWallSize = sizeList.get(sizeList.size() / 2);
medianWallSize = medianWallSize * 2 + 1;
int maxWallSize = (int) (medianWallSize * 1.5);
// Mask out erroneously dilated pixels
double avgx = 0;
double avgy = 0;
double count = 0;
for (int x = startX; x < stopX; x++) {
for (int y = startY; y < stopY; y++) {
int i = x + y * width;
if (currentWall[i] != 0 && pixels[i] != (short) 255 // not a foreground pixel
|| lumen.getDistanceMap()[i] > maxWallSize) { // outside max wallsize
currentWall[i] = 0;
}
// remove lumenmark from data
if (lumen.getLumenMask()[i] != 0) {
count++;
avgx += x;
avgy += y;
pixels[i] = 0;
}
if (currentWall[i] != 0) {
lumen.getLumenMask()[i] = (short) 255;
}
}
}
sout[(int) (avgx / count) + (int) (avgy / count) * width] = lumen.getId();
}
skeletonLog.putSlice(sout, z);
for (int i = 0; i < pixels.length; i++) {
int min = Integer.MAX_VALUE;
short minId = 0;
if (pixels[i] != 0) {
for (int j = 0; j < lumens.size(); j++) {