package Algorithms;
import Utils.*;
import java.util.Arrays;
public class AnalyzeSkeleton2
{
public static final int NONE = 0;
public static final int END_POINT = 1;
public static final int JUNCTION = 2;
public static final int SLAB = 3;
public static final int MAX_BREADTH = 8;
public static final int JUNC_VISIT = 16;
public static final int SKEL_VISIT = 24;
public static class Result
{
public IntVector[] pointVectors;
public IntVector endPoints;
public IntVector junctionVoxels;
public IntVector singleJunctions;
public IntVector slabVoxels;
public IntVector startingSlabVoxels;
public IntVector toRevisit;
public IntVector isolatedJunctions;
public IntVector removedJunctions;
public IntVector slabList;
public IntVector edgesTrees;
public IntVector edgesVerts;
public IntVector edgesPoints;
public DoubleVector edgesLengths;
public int treeCount;
//public IntVector triplePointCounts;
//public IntVector quadruplePointCounts;
public DoubleVector totalBranchLengths;
public DoubleVector maximumBranchLengths;
public IntVector numberOfBranches;
public IntVector numberOfSlabs;
private static IntVector resetIntVector(IntVector vec)
{
if (vec == null)
return new IntVector();
vec.size = 0;
return vec;
}
private static DoubleVector resetDoubleVector(DoubleVector vec)
{
if (vec == null)
return new DoubleVector();
vec.size = 0;
return vec;
}
public void reset()
{
if (pointVectors == null)
pointVectors = new IntVector[3];
endPoints = resetIntVector(endPoints);
junctionVoxels = resetIntVector(junctionVoxels);
singleJunctions = resetIntVector(singleJunctions);
slabVoxels = resetIntVector(slabVoxels);
startingSlabVoxels = resetIntVector(startingSlabVoxels);
toRevisit = resetIntVector(toRevisit);
isolatedJunctions = resetIntVector(isolatedJunctions);
removedJunctions = resetIntVector(removedJunctions);
slabList = resetIntVector(slabList);
edgesTrees = resetIntVector(edgesTrees);
edgesVerts = resetIntVector(edgesVerts);
edgesPoints = resetIntVector(edgesPoints);
edgesLengths = resetDoubleVector(edgesLengths);
treeCount = 0;
//triplePointCounts = resetIntVector(triplePointCounts);
//quadruplePointCounts = resetIntVector(quadruplePointCounts);
totalBranchLengths = resetDoubleVector(totalBranchLengths);
maximumBranchLengths = resetDoubleVector(maximumBranchLengths);
numberOfBranches = resetIntVector(numberOfBranches);
numberOfSlabs = resetIntVector(numberOfSlabs);
}
}
private static class Scratch
{
public int width;
public int height;
public int breadth;
public int[] imageInfo;
public int[] junctionMap2d;
public int[] markedImages; // volume
public int[] endPointVertexMap; // volume
public int[] junctionVertexMap; // volume
public void acquire(int width, int height, int breadth)
{
this.width = width;
this.height = height;
this.breadth = breadth;
int area = width * height;
this.imageInfo = IntBufferPool.acquireAsIs(area);
this.junctionMap2d = IntBufferPool.acquireZeroed(area);
this.markedImages = IntBufferPool.acquireZeroed(area * breadth);
this.endPointVertexMap = IntBufferPool.acquireZeroed(area * breadth);
this.junctionVertexMap = IntBufferPool.acquireZeroed(area * breadth);
}
public void release()
{
this.imageInfo = IntBufferPool.release(imageInfo);
this.junctionMap2d = IntBufferPool.release(junctionMap2d);
this.markedImages = IntBufferPool.release(markedImages);
this.endPointVertexMap = IntBufferPool.release(endPointVertexMap);
this.junctionVertexMap = IntBufferPool.release(junctionVertexMap);
}
}
// returns a 26-bit number containing each neighbor within a 3x3x3 vicinity (-1 to exclude the point itself)
public static int getBooleanNeighborBits(byte[] planes, int width, int height, int breadth, int x, int y, int z)
{
int bits = 0;
for (int i = 0; i < 27; i++) {
if (i == 13)
continue;
int xx = x + (i % 3) - 1;
int yy = y + ((i / 3) % 3) - 1;
int zz = z + (i / 9) - 1;
int p = 0;
if (xx >= 0 && xx < width && yy >= 0 && yy < height && zz >= 0 && zz < breadth)
p = (planes[xx + width * yy] >> zz) & 1;
bits = (bits << 1) | p;
}
return bits;
}
public static void analyze(
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth
) {
result.reset();
breadth = Math.min(Math.max(1, breadth), MAX_BREADTH);
Scratch scratch = new Scratch();
scratch.acquire(width, height, breadth);
tagImages(scratch, result, skeletonImages, width, height, breadth);
int nTrees = markTreesOnImages(scratch, result, skeletonImages, width, height, breadth);
if (nTrees <= 0)
nTrees = 1;
result.treeCount = nTrees;
//resizeAndZeroIntVector(result.triplePointCounts, nTrees);
//resizeAndZeroIntVector(result.quadruplePointCounts, nTrees);
groupJunctions(scratch, result, skeletonImages, width, height, breadth, nTrees);
resizeAndZeroDoubleVector(result.totalBranchLengths, nTrees);
resizeAndZeroDoubleVector(result.maximumBranchLengths, nTrees);
resizeAndZeroIntVector(result.numberOfBranches, nTrees);
resizeAndZeroIntVector(result.numberOfSlabs, nTrees);
buildSkeletonGraphs(scratch, result, skeletonImages, width, height, breadth, nTrees);
isolateDominantJunctions(scratch, result, width, height);
scratch.release();
}
static void tagImages(
Scratch scratch,
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth
) {
result.pointVectors[END_POINT-1] = result.endPoints;
result.pointVectors[JUNCTION-1] = result.junctionVoxels;
result.pointVectors[SLAB-1] = result.slabVoxels;
int area = width * height;
for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
int idx = x + width * y;
int value = 0;
for (int z = 0; z < breadth; ++z) {
int type = NONE;
if (((skeletonImages[idx] >>> z) & 1) != 0) {
int numOfNeighbors = Integer.bitCount(
getBooleanNeighborBits(skeletonImages, width, height, breadth, x, y, z)
);
if (numOfNeighbors < 2) {
type = END_POINT;
int pos = result.pointVectors[type-1].addThree(x, y, z);
scratch.endPointVertexMap[idx + area * z] = pos;
}
else if (numOfNeighbors > 2) {
type = JUNCTION;
int pos = result.pointVectors[type-1].addThree(x, y, z);
scratch.junctionMap2d[idx] = pos + 1;
}
else {
type = SLAB;
result.pointVectors[type-1].addThree(x, y, z);
}
}
value = (value << 2) | type;
}
scratch.imageInfo[idx] = value;
}
}
}
static int markTreesOnImages(
Scratch scratch,
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth
) {
final int area = width * height;
int nTrees = 0;
for (int type = END_POINT; type <= SLAB; type++) {
int n = result.pointVectors[type-1].size;
for (int i = 0; i < n; i += 3) {
int x = result.pointVectors[type-1].buf[i];
int y = result.pointVectors[type-1].buf[i+1];
int z = result.pointVectors[type-1].buf[i+2];
if (scratch.markedImages[x + width * y + area * z] != 0)
continue;
if (type == SLAB)
result.startingSlabVoxels.addThree(x, y, z);
int numOfVoxels = 0;
scratch.markedImages[x + width * y + area * z] = nTrees + 1;
result.toRevisit.size = 0;
if (((scratch.imageInfo[x + width * y] >>> (z << 1)) & 3) == JUNCTION)
result.toRevisit.addThree(x, y, z);
int revisitIdx = 0;
boolean wasRevisit = false;
boolean didFindPoint;
do {
didFindPoint = false;
for (int j = 0; j < 27; j++) {
if (j == 13)
continue;
int xx = x + (j / 9) - 1;
int yy = y + ((j / 3) % 3) - 1;
int zz = z + (j % 3) - 1;
if (xx < 0 || xx >= width || yy < 0 || yy >= height || zz < 0 || zz >= breadth)
continue;
if (
((skeletonImages[xx + width * yy] >>> zz) & 1) != 0 &&
scratch.markedImages[xx + width * yy + area * zz] == 0
) {
x = xx;
y = yy;
z = zz;
didFindPoint = true;
break;
}
}
if (didFindPoint) {
++numOfVoxels;
scratch.markedImages[x + width * y + area * z] = nTrees + 1;
if (((scratch.imageInfo[x + width * y] >>> (z << 1)) & 3) == JUNCTION)
result.toRevisit.addThree(x, y, z);
}
else {
if (wasRevisit)
revisitIdx += 3;
if (revisitIdx < result.toRevisit.size) {
x = result.toRevisit.buf[revisitIdx];
y = result.toRevisit.buf[revisitIdx+1];
z = result.toRevisit.buf[revisitIdx+2];
}
}
wasRevisit = !didFindPoint;
}
while (didFindPoint || revisitIdx < result.toRevisit.size);
if (type == END_POINT || numOfVoxels != 0)
nTrees++;
}
}
return nTrees;
}
static void groupJunctions(
Scratch scratch,
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth,
int nTrees
) {
final int area = width * height;
int vertexCount = 0;
for (int i = 0; i < result.junctionVoxels.size; i += 3) {
int x = result.junctionVoxels.buf[i];
int y = result.junctionVoxels.buf[i+1];
int z = result.junctionVoxels.buf[i+2];
if (((scratch.imageInfo[x + width * y] >>> (JUNC_VISIT + z)) & 1) != 0)
continue;
scratch.imageInfo[x + width * y] |= 1 << (JUNC_VISIT + z);
int treeIdx = scratch.markedImages[x + width * y + area * z] - 1;
vertexCount++;
scratch.junctionVertexMap[x + width * y + area * z] = vertexCount;
result.singleJunctions.addThree(x, y, z);
result.toRevisit.size = 0;
result.toRevisit.addThree(x, y, z);
int revisitIdx = 0;
boolean wasRevisit = false;
int pointFound;
do {
pointFound = -1;
int nBranches = 0;
for (int j = 0; j < 27; j++) {
if (j == 13)
continue;
int xx = x + (j / 9) - 1;
int yy = y + ((j / 3) % 3) - 1;
int zz = z + (j % 3) - 1;
if (xx < 0 || xx >= width || yy < 0 || yy >= height || zz < 0 || zz >= breadth)
continue;
int type = (scratch.imageInfo[xx + width * yy] >>> (zz << 1)) & 3;
if (
pointFound == -1 &&
((skeletonImages[xx + width * yy] >>> zz) & 1) != 0 &&
((scratch.imageInfo[xx + width * yy] >>> (JUNC_VISIT + zz)) & 1) == 0 &&
type == JUNCTION
) {
pointFound = j;
}
if (type == END_POINT || type == SLAB)
nBranches++;
}
if (pointFound >= 0) {
x += (pointFound / 9) - 1;
y += ((pointFound / 3) % 3) - 1;
z += (pointFound % 3) - 1;
scratch.imageInfo[x + width * y] |= 1 << (JUNC_VISIT + z);
scratch.junctionVertexMap[x + width * y + area * z] = vertexCount;
result.singleJunctions.addThree(x, y, z);
result.toRevisit.addThree(x, y, z);
/*
if (nBranches == 3)
result.triplePointCounts.buf[treeIdx]++;
else if (nBranches == 4)
result.quadruplePointCounts.buf[treeIdx]++;
*/
}
else {
if (wasRevisit)
revisitIdx += 3;
if (revisitIdx < result.toRevisit.size) {
x = result.toRevisit.buf[revisitIdx];
y = result.toRevisit.buf[revisitIdx+1];
z = result.toRevisit.buf[revisitIdx+2];
}
}
wasRevisit = pointFound == -1;
}
while (pointFound >= 0 || revisitIdx < result.toRevisit.size);
}
}
static void buildSkeletonGraphs(
Scratch scratch,
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth,
int nTrees
) {
int area = width * height;
for (int i = 0; i < result.endPoints.size; i += 3) {
int x = result.endPoints.buf[i];
int y = result.endPoints.buf[i+1];
int z = result.endPoints.buf[i+2];
int t = scratch.markedImages[x + width * y + area * z] - 1;
if (((scratch.imageInfo[x + width * y] >>> (SKEL_VISIT + z)) & 1) != 0)
continue;
scratch.imageInfo[x + width * y] |= 1 << (SKEL_VISIT + z);
int slabListIdx = result.slabList.size;
visitBranch(
scratch, result, skeletonImages, width, height, breadth, t,
END_POINT, i, slabListIdx, 0.0, x, y, z
);
}
for (int i = 0; i < result.singleJunctions.size; i += 3) {
int x = result.singleJunctions.buf[i];
int y = result.singleJunctions.buf[i+1];
int z = result.singleJunctions.buf[i+2];
int t = scratch.markedImages[x + width * y + area * z] - 1;
// no check
scratch.imageInfo[x + width * y] |= 1 << (SKEL_VISIT + z);
int vertexIdx = scratch.junctionVertexMap[x + width * y + area * z] - 1;
for (int j = 0; j < 27; j++) {
if (j == 13)
continue;
int xx = x + (j / 9) - 1;
int yy = y + ((j / 3) % 3) - 1;
int zz = z + (j % 3) - 1;
if (xx < 0 || xx >= width || yy < 0 || yy >= height || zz < 0 || zz >= breadth)
continue;
if (
((skeletonImages[xx + width * yy] >>> zz) & 1) != 0 &&
((scratch.imageInfo[xx + width * yy] >>> (SKEL_VISIT + zz)) & 1) == 0
) {
//didFindPoint = true;
if (((scratch.imageInfo[xx + width * yy] >>> (zz << 1)) & 3) == JUNCTION) {
scratch.imageInfo[xx + width * yy] |= 1 << (SKEL_VISIT + zz);
continue;
}
double initialLength = calculateDistance(x, y, z, xx, yy, zz);
int slabListIdx = result.slabList.addThree(xx, yy, zz);
visitBranch(
scratch, result, skeletonImages, width, height, breadth, t,
JUNCTION, vertexIdx, slabListIdx, initialLength, xx, yy, zz
);
}
}
}
int startingSlabStart = 0;
while (startingSlabStart + 2 < result.startingSlabVoxels.size) {
int s = startingSlabStart;
int x = result.startingSlabVoxels.buf[s];
int y = result.startingSlabVoxels.buf[s+1];
int z = result.startingSlabVoxels.buf[s+2];
int t = scratch.markedImages[x + width * y + area * z] - 1;
boolean isSingle = true;
startingSlabStart += 3;
while (startingSlabStart + 2 < result.startingSlabVoxels.size) {
int xx = result.startingSlabVoxels.buf[startingSlabStart];
int yy = result.startingSlabVoxels.buf[startingSlabStart+1];
int zz = result.startingSlabVoxels.buf[startingSlabStart+2];
if (scratch.markedImages[xx + width * yy + area * zz] != t + 1)
break;
isSingle = false;
startingSlabStart += 3;
}
if (isSingle) {
//System.out.println("slab visit on tree " + (t+1) + " / " + nTrees);
result.numberOfSlabs.buf[t]++;
int slabListIdx = result.slabList.addThree(x, y, z);
visitBranch(
scratch, result, skeletonImages, width, height, breadth, t,
SLAB, s, slabListIdx, 0.0, x, y, z
);
}
}
}
static void visitBranch(
Scratch scratch,
Result result,
byte[] skeletonImages,
int width,
int height,
int breadth,
int iTree,
int mode,
int initialVertIdx,
int initialSlabListIdx,
double initialLength,
int xStart,
int yStart,
int zStart
) {
final int area = width * height;
double length = initialLength;
int type = NONE;
int finalVertIdx = initialVertIdx;
int x = xStart;
int y = yStart;
int z = zStart;
scratch.imageInfo[x + width * y] |= 1 << (SKEL_VISIT + z);
boolean didFindSlabPoint;
do {
didFindSlabPoint = false;
for (int j = 0; j < 27; j++) {
if (j == 13)
continue;
int xx = x + (j / 9) - 1;
int yy = y + ((j / 3) % 3) - 1;
int zz = z + (j % 3) - 1;
if (xx < 0 || xx >= width || yy < 0 || yy >= height || zz < 0 || zz >= breadth)
continue;
if (
((skeletonImages[xx + width * yy] >>> zz) & 1) != 0 &&
((scratch.imageInfo[xx + width * yy] >>> (SKEL_VISIT + zz)) & 1) == 0
) {
length += calculateDistance(x, y, z, xx, yy, zz);
scratch.imageInfo[xx + width * yy] |= 1 << (SKEL_VISIT + zz);
type = (scratch.imageInfo[xx + width * yy] >>> (zz << 1)) & 3;
x = xx;
y = yy;
z = zz;
switch (type) {
case END_POINT:
finalVertIdx = scratch.endPointVertexMap[x + width * y + area * z] - 1;
break;
case JUNCTION:
finalVertIdx = scratch.junctionVertexMap[x + width * y + area * z] - 1;
//this.auxFinalVertex = this.findPointVertex(this.junctionVertex[iTree], nextPoint);
break;
case SLAB:
result.numberOfSlabs.buf[iTree]++;
result.slabList.addThree(x, y, z);
didFindSlabPoint = true;
break;
}
break;
}
}
}
while (didFindSlabPoint);
double lengthBeforeSlabs = length;
int modeEnd = mode;
if (mode == SLAB) {
finalVertIdx = initialVertIdx;
}
else if (type == SLAB) {
finalVertIdx = initialVertIdx;
int xExclude = -1;
int yExclude = -1;
int zExclude = -1;
int vertExclude = -1;
if (mode == END_POINT) {
xExclude = xStart;
yExclude = yStart;
zExclude = zStart;
}
else {
vertExclude = initialVertIdx;
}
for (int j = 0; j < 27; j++) {
if (j == 13)
continue;
int xx = x + (j / 9) - 1;
int yy = y + ((j / 3) % 3) - 1;
int zz = z + (j % 3) - 1;
if (xx < 0 || xx >= width || yy < 0 || yy >= height || zz < 0 || zz >= breadth)
continue;
int info = scratch.imageInfo[xx + width * yy];
int vert = scratch.junctionVertexMap[xx + width * yy + area * zz];
if (
((skeletonImages[xx + width * yy] >>> zz) & 1) != 0 &&
((info >>> (SKEL_VISIT + zz)) & 1) != 0 &&
((info >>> (zz << 1)) & 3) == JUNCTION &&
vert != vertExclude &&
!(xx == xExclude && yy == yExclude && zz == zExclude)
) {
finalVertIdx = vert;
modeEnd = JUNCTION;
length += calculateDistance(x, y, z, xx, yy, zz);
break;
}
}
}
result.totalBranchLengths.buf[iTree] += mode == JUNCTION ? lengthBeforeSlabs : length;
if (length > result.maximumBranchLengths.buf[iTree])
result.maximumBranchLengths.buf[iTree] = length;
result.numberOfBranches.buf[iTree]++;
int modeSides = (mode << 2) | modeEnd;
int nSlabs = result.slabList.size - initialSlabListIdx;
result.edgesTrees.add(iTree);
result.edgesVerts.addThree(modeSides, initialVertIdx, finalVertIdx);
result.edgesPoints.addTwo(initialSlabListIdx, nSlabs);
result.edgesLengths.add(length);
}
static void isolateDominantJunctions(Scratch scratch, Result result, int width, int height)
{
for (int i = 0; i < result.junctionVoxels.size; i += 3) {
int x = result.junctionVoxels.buf[i];
int y = result.junctionVoxels.buf[i+1];
int pos = scratch.junctionMap2d[x + width * y];
if (pos <= 0)
continue;
result.isolatedJunctions.add(pos - 1);
scratch.junctionMap2d[x + width * y] = 0;
for (int j = 0; j < 9; j++) {
if (j == 4)
continue;
int xx = x + (j % 3) - 1;
int yy = y + (j / 3) - 1;
int p = scratch.junctionMap2d[xx + width * yy];
if (p > 0) {
result.removedJunctions.add(p - 1);
scratch.junctionMap2d[xx + width * yy] = 0;
}
}
}
}
static double calculateDistance(int x1, int y1, int z1, int x2, int y2, int z2)
{
double dx = (double)(x2 - x1);
double dy = (double)(y2 - y1);
double dz = (double)(z2 - z1);
return Math.sqrt((dx * dx) + (dy * dy) + (dz * dz));
}
static void resizeAndZeroDoubleVector(DoubleVector vec, int newSize)
{
vec.resize(newSize);
if (newSize > 0)
Arrays.fill(vec.buf, 0, newSize, 0.0);
}
static void resizeAndZeroIntVector(IntVector vec, int newSize)
{
vec.resize(newSize);
if (newSize > 0)
Arrays.fill(vec.buf, 0, newSize, 0);
}
}