Example 1 with FastMarchingMethodHA
use of cbit.vcell.render.FastMarchingMethodHA in project vcell by virtualcell.
the class DistanceMapGenerator method computeDistanceMaps.
public static SubvolumeSignedDistanceMap[] computeDistanceMaps(Geometry geometry, VCImage subvolumeHandleImage, boolean bCellCentered, boolean insideOnly) throws ImageException {
double[] samplesX = new double[subvolumeHandleImage.getNumX()];
double[] samplesY = new double[subvolumeHandleImage.getNumY()];
double[] samplesZ = new double[subvolumeHandleImage.getNumZ()];
ISize sampleSize = new ISize(subvolumeHandleImage.getNumX(), subvolumeHandleImage.getNumY(), subvolumeHandleImage.getNumZ());
byte[] pixels = subvolumeHandleImage.getPixels();
boolean[] ignoreMask = new boolean[sampleSize.getXYZ()];
Origin origin = geometry.getOrigin();
Extent extent = geometry.getExtent();
RayCaster.sampleXYZCoordinates(sampleSize, origin, extent, samplesX, samplesY, samplesZ, bCellCentered);
ArrayList<SubvolumeSignedDistanceMap> distanceMaps = new ArrayList<SubvolumeSignedDistanceMap>();
int count = 0;
for (SubVolume subVolume : geometry.getGeometrySpec().getSubVolumes()) {
//
// find surfaces that bound the current SubVolume
//
ArrayList<Surface> surfaces = new ArrayList<Surface>();
for (GeometricRegion geometricRegion : geometry.getGeometrySurfaceDescription().getGeometricRegions()) {
if (geometricRegion instanceof SurfaceGeometricRegion) {
SurfaceGeometricRegion surfaceGeometricRegion = (SurfaceGeometricRegion) geometricRegion;
for (GeometricRegion adjacentRegion : ((SurfaceGeometricRegion) geometricRegion).getAdjacentGeometricRegions()) {
if (adjacentRegion instanceof VolumeGeometricRegion && ((VolumeGeometricRegion) adjacentRegion).getSubVolume() == subVolume) {
surfaces.add(geometry.getGeometrySurfaceDescription().getSurfaceCollection().getSurface(surfaceGeometricRegion));
}
}
}
}
// find unsigned distances in a narrow band for surfaces that bound this SubVolume (expensive)
// values outside the band are assumed to be initialized to MAX_NUMBER
long t1 = System.currentTimeMillis();
double[] distanceMap = localUnsignedDistanceMap(surfaces, samplesX, samplesY, samplesZ, 3);
long t2 = System.currentTimeMillis();
System.out.println(" Distance to triangle: " + (int) ((t2 - t1) / 1000) + " sec.");
// extend signed distance map using fast marching method from narrow band to all points.
// will do it in 2 steps, positive growth first towards inside, then change the sign of the whole
// distance map, then positive growth towards the exterior
// this way, the interior distances will end negative and the exterior distances positive
// 2 step growth saves memory and reduces the number of elements present at any given time in the binary
// heap (binary heap manipulation is the most time consuming factor and it depends on the # of elements)
Arrays.fill(ignoreMask, true);
int subvolumeHandle = subVolume.getHandle();
for (int i = 0; i < ignoreMask.length; i++) {
if (pixels[i] == subvolumeHandle) {
// inside
ignoreMask[i] = false;
} else {
// outside
if (distanceMap[i] < MAX_NUMBER) {
// make negative the part of narrow band which is outside
distanceMap[i] = -distanceMap[i];
}
}
}
// // step 1, we compute distances for the points "inside"
// // the points outside are cold (we don't compute their distances this step)
double deltaX = samplesX[1] - samplesX[0];
double deltaY = samplesY[1] - samplesY[0];
double deltaZ = samplesZ[1] - samplesZ[0];
FastMarchingMethodHA fmm = new FastMarchingMethodHA(samplesX.length, samplesY.length, samplesZ.length, deltaX, deltaY, deltaZ, distanceMap, ignoreMask);
fmm.march();
if (!insideOnly) {
// sign change of the half-completed distance map, the "interior" will become negative as it should be
for (int i = 0; i < distanceMap.length; i++) {
if (distanceMap[i] < MAX_NUMBER) {
distanceMap[i] = -distanceMap[i];
}
}
// step 2, we compute distances for the points "outside"
// no cold points (points we don't care about) this time, they are already frozen
fmm = new FastMarchingMethodHA(samplesX.length, samplesY.length, samplesZ.length, deltaX, deltaY, deltaZ, distanceMap, null);
fmm.march();
} else {
// sign change of the half-completed distance map, the "interior" will become negative as it should be
for (int i = 0; i < distanceMap.length; i++) {
if (distanceMap[i] < MAX_NUMBER) {
if (pixels[i] != subvolumeHandle) {
// need to filter out the part of the narrow band which is not inside
distanceMap[i] = MAX_NUMBER;
} else {
distanceMap[i] = -distanceMap[i];
}
}
}
}
// try { // save some points in a VisIt compatible format
// int xm = samplesX.length;
// int ym = samplesY.length;
// BufferedWriter out = new BufferedWriter(new FileWriter("c:\\TEMP\\2D_circle" + count + ".3D"));
// out.write("x y z value\n");
//
// for(int j=0; j<distanceMap.length; j++) {
// int x = getX(j,xm,ym);
// int y = getY(j,xm,ym);
// int z = getZ(j,xm,ym);
// if(distanceMap[j] < MAX_NUMBER) {
// if((j%2 == 0 || j%3 == 0) && (distanceMap[j] <= -2)) {
// out.write(x + " " + y + " " + z + " " + (int)(distanceMap[j]*10) + "\n");
// } else if((j%17 == 0 || j%23 == 0) && (distanceMap[j] <= 0.5) && (distanceMap[j] > -2)) {
// out.write(x + " " + y + " " + z + " " + (int)(distanceMap[j]*10) + "\n");
// } else if((j%31 == 0 || j%41 == 0) && (distanceMap[j] > 0.5)) {
// out.write(x + " " + y + " " + z + " " + (int)(distanceMap[j]*10) + "\n");
// }
// }
// if(distanceMap[j] < MAX_NUMBER) {
// if(j%2 == 0) {
// out.write(x + " " + y + " " + z + " " + (int)(distanceMap[j]*100) + "\n");
// }
// }
// if(x==50 && y==50 && z==25) { // on the surface
// out.write(x + " " + y + " " + z + " " + (int)(distanceMap[j]*100) + "\n");
// }
// if(x==0 && y==0 && z==0) {
// out.write(x + " " + y + " " + z + " " + 0 + "\n");
// }
// if(x==100 && y==100 && z==100) {
// out.write(x + " " + y + " " + z + " " + 0 + "\n");
// }
//
// }
// out.close();
// } catch (IOException e) {
// }
SubvolumeSignedDistanceMap subvolumeSignedDistanceMap = new SubvolumeSignedDistanceMap(subVolume, samplesX, samplesY, samplesZ, distanceMap);
distanceMaps.add(subvolumeSignedDistanceMap);
count++;
}
return distanceMaps.toArray(new SubvolumeSignedDistanceMap[distanceMaps.size()]);
}