Examples with RegionImage cbit.vcell.geometry.RegionImage used on opensource projects

Example 21 with RegionImage

use of cbit.vcell.geometry.RegionImage in project vcell by virtualcell.

the class GeometrySurfaceDescription method updateAll.

/**
 * Insert the method's description here.
 * Creation date: (5/26/2004 10:19:59 AM)
 */
public void updateAll() throws GeometryException, ImageException, ExpressionException {
    // 
    // updates if necessary: regionImage, surfaceCollection and resolvedLocations[]
    // 
    // assumes that volumeSampleSize and filterCutoffFrequency are already specified
    // 
    // 
    // 
    // make new RegionImage if necessary missing or wrong size
    // 
    VCellThreadChecker.checkCpuIntensiveInvocation();
    boolean bChanged = false;
    RegionImage updatedRegionImage = GeometrySurfaceUtils.getUpdatedRegionImage(this);
    if (updatedRegionImage != getRegionImage()) {
        // getUpdatedRegionImage returns same image if not changed
        getRegionImage0().setValue(updatedRegionImage);
        bChanged = true;
    }
    // 
    if (getSurfaceCollection() == null || bChanged) {
        setSurfaceCollection(updatedRegionImage.getSurfacecollection());
        bChanged = true;
    }
    // 
    if (getGeometricRegions() == null || bChanged) {
        try {
            setGeometricRegions(GeometrySurfaceUtils.getUpdatedGeometricRegions(this, getRegionImage(), getSurfaceCollection()));
            bChanged = true;
        } catch (java.beans.PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new GeometryException("unexpected exception while generating regions: " + e.getMessage());
        }
    }
    // 
    if (getSurfaceClasses() == null || bChanged) {
        refreshSurfaceClasses();
        bChanged = true;
    }
}

Example 22 with RegionImage

use of cbit.vcell.geometry.RegionImage in project vcell by virtualcell.

the class GeometrySurfaceUtils method getUpdatedGeometricRegions.

/**
 * Insert the method's description here.
 * Creation date: (6/28/2004 3:09:13 PM)
 * @return cbit.vcell.geometry.surface.GeometricRegion[]
 * @param geoSurfaceDescription cbit.vcell.geometry.surface.GeometrySurfaceDescription
 * @param surfaceCollection cbit.vcell.geometry.surface.SurfaceCollection
 */
public static GeometricRegion[] getUpdatedGeometricRegions(GeometrySurfaceDescription geoSurfaceDescription, cbit.vcell.geometry.RegionImage regionImage, SurfaceCollection surfaceCollection) {
    // 
    // parse regionImage into ResolvedVolumeLocations
    // 
    VCellThreadChecker.checkCpuIntensiveInvocation();
    double sizeOfPixel = 0;
    VCUnitDefinition volumeUnit = null;
    VCUnitDefinition surfaceUnit = null;
    GeometrySpec geometrySpec = geoSurfaceDescription.getGeometry().getGeometrySpec();
    GeometryUnitSystem geometryUnitSystem = geoSurfaceDescription.getGeometry().getUnitSystem();
    switch(geometrySpec.getDimension()) {
        case 1:
            {
                sizeOfPixel = geometrySpec.getExtent().getX() / (regionImage.getNumX() - 1);
                // sizeOfPixel /= 9.0;  // to account for the padding from 1D to 3D
                volumeUnit = geometryUnitSystem.getLengthUnit();
                surfaceUnit = geometryUnitSystem.getInstance_DIMENSIONLESS();
                break;
            }
        case 2:
            {
                sizeOfPixel = geometrySpec.getExtent().getX() / (regionImage.getNumX() - 1) * geometrySpec.getExtent().getY() / (regionImage.getNumY() - 1);
                // sizeOfPixel /= 3.0;  // to account for the padding from 2D to 3D
                volumeUnit = geometryUnitSystem.getAreaUnit();
                surfaceUnit = geometryUnitSystem.getLengthUnit();
                break;
            }
        case 3:
            {
                sizeOfPixel = geometrySpec.getExtent().getX() / (regionImage.getNumX() - 1) * geometrySpec.getExtent().getY() / (regionImage.getNumY() - 1) * geometrySpec.getExtent().getZ() / (regionImage.getNumZ() - 1);
                volumeUnit = geometryUnitSystem.getVolumeUnit();
                surfaceUnit = geometryUnitSystem.getAreaUnit();
                break;
            }
    }
    int numX = regionImage.getNumX();
    int numY = regionImage.getNumY();
    int numZ = regionImage.getNumZ();
    int numXY = numX * numY;
    java.util.Vector<GeometricRegion> regionList = new java.util.Vector<GeometricRegion>();
    cbit.vcell.geometry.RegionImage.RegionInfo[] regionInfos = regionImage.getRegionInfos();
    for (int i = 0; i < regionInfos.length; i++) {
        cbit.vcell.geometry.RegionImage.RegionInfo regionInfo = regionInfos[i];
        lg.info(regionInfo);
        cbit.vcell.geometry.SubVolume subVolume = geometrySpec.getSubVolume(regionInfo.getPixelValue());
        String name = subVolume.getName() + regionInfo.getRegionIndex();
        int numPixels = regionInfo.getNumPixels();
        double size = numPixels * sizeOfPixel;
        // 
        switch(geometrySpec.getDimension()) {
            case 1:
                {
                    int numHalvesToRemove = 0;
                    // -x pixel
                    if (regionInfo.isIndexInRegion(0)) {
                        // (regionImage.getRegionIndex(0,0,0) == regionIndex)
                        numHalvesToRemove++;
                    }
                    // +x pixel
                    if (regionInfo.isIndexInRegion(numX - 1)) {
                        // (regionImage.getRegionIndex(regionImageNumX-1,0,0) == regionIndex){
                        numHalvesToRemove++;
                    }
                    // *9; // 1D is padded by 3 in y and 3 in z (hence factor of 9)
                    size -= sizeOfPixel * 0.5 * numHalvesToRemove;
                    break;
                }
            case 2:
                {
                    int numQuadrantsToRemove = 0;
                    int yOffset = 0;
                    for (int yIndex = 0; yIndex < numY; yIndex++) {
                        if (regionInfo.isIndexInRegion(yOffset)) {
                            // (regionImage.getRegionIndex(0,yIndex,0) == regionIndex){
                            if (yIndex == 0 || yIndex == numY - 1) {
                                // corner, remove 3 quadrants
                                numQuadrantsToRemove += 3;
                            } else {
                                // side, remove 2 quadrants
                                numQuadrantsToRemove += 2;
                            }
                        }
                        // +x side (including attached corners)
                        if (regionInfo.isIndexInRegion(numX - 1 + yOffset)) {
                            // (regionImage.getRegionIndex(numX-1,yIndex,0) == regionIndex){
                            if (yIndex == 0 || yIndex == numY - 1) {
                                // corner, remove 3 quadrants
                                numQuadrantsToRemove += 3;
                            } else {
                                // side, remove 2 quadrants
                                numQuadrantsToRemove += 2;
                            }
                        }
                        yOffset += numX;
                    }
                    int yOffsetLastLine = (numY - 1) * numX;
                    for (int xIndex = 1; xIndex < numX - 1; xIndex++) {
                        // -y side (excluding corners)
                        if (regionInfo.isIndexInRegion(xIndex)) {
                            // (regionImage.getRegionIndex(xIndex,0,0) == regionIndex){
                            // side, remove 2 quadrants
                            numQuadrantsToRemove += 2;
                        }
                        // +x side (excluding corners)
                        if (regionInfo.isIndexInRegion(xIndex + yOffsetLastLine)) {
                            // (regionImage.getRegionIndex(xIndex,numY-1,0) == regionIndex){
                            // side, remove 2 quadrants
                            numQuadrantsToRemove += 2;
                        }
                    }
                    // *3; // 2D is padded by 3 in z (hence the factor of 3).
                    size -= sizeOfPixel * 0.25 * numQuadrantsToRemove;
                    break;
                }
            case 3:
                {
                    int numOctantsToRemove = 0;
                    for (int zIndex = 0; zIndex < numZ; zIndex++) {
                        for (int yIndex = 0; yIndex < numY; yIndex++) {
                            // -x side (including attached edges and corners)
                            // already on face of boundary (removing half)
                            int totalOctants = 4;
                            // lg.info("-x side including edges and corners");
                            if (regionInfo.isIndexInRegion(yIndex * numX + zIndex * numXY)) {
                                // (regionImage.getRegionIndex(0,yIndex,zIndex) == regionIndex){
                                if (yIndex == 0 || yIndex == numY - 1) {
                                    totalOctants /= 2;
                                }
                                if (zIndex == 0 || zIndex == numZ - 1) {
                                    totalOctants /= 2;
                                }
                                numOctantsToRemove += (8 - totalOctants);
                            }
                            // lg.info("+x side including edges and corners");
                            // +x side (including attached edges and corners)
                            // already on face of boundary (removing half)
                            totalOctants = 4;
                            if (regionInfo.isIndexInRegion(numX - 1 + yIndex * numX + zIndex * numXY)) {
                                // (regionImage.getRegionIndex(numX-1,yIndex,zIndex) == regionIndex){
                                if (yIndex == 0 || yIndex == numY - 1) {
                                    totalOctants /= 2;
                                }
                                if (zIndex == 0 || zIndex == numZ - 1) {
                                    totalOctants /= 2;
                                }
                                numOctantsToRemove += (8 - totalOctants);
                            }
                        }
                        for (int xIndex = 1; xIndex < numX - 1; xIndex++) {
                            // -y side (including attached edges along x axis, excluding corners)
                            // lg.info("-y side including edges");
                            // already on face of boundary (removing half)
                            int totalOctants = 4;
                            if (regionInfo.isIndexInRegion(xIndex + zIndex * numXY)) {
                                // (regionImage.getRegionIndex(xIndex,0,zIndex) == regionIndex){
                                if (zIndex == 0 || zIndex == numZ - 1) {
                                    totalOctants /= 2;
                                }
                                numOctantsToRemove += (8 - totalOctants);
                            }
                            // lg.info("+y side including edges");
                            // +y side (including attached edges along x axis, excluding corners)
                            // already on face of boundary (removing half)
                            totalOctants = 4;
                            if (regionInfo.isIndexInRegion(xIndex + (numY - 1) * numX + zIndex * numXY)) {
                                // (regionImage.getRegionIndex(xIndex,numY-1,zIndex) == regionIndex){
                                if (zIndex == 0 || zIndex == numZ - 1) {
                                    totalOctants /= 2;
                                }
                                numOctantsToRemove += (8 - totalOctants);
                            }
                        }
                    }
                    for (int yIndex = 1; yIndex < numY - 1; yIndex++) {
                        for (int xIndex = 1; xIndex < numX - 1; xIndex++) {
                            // already on face of boundary (removing half)
                            int totalOctants = 4;
                            // lg.info("-z side including nothing");
                            if (regionInfo.isIndexInRegion(xIndex + yIndex * numX)) {
                                // (regionImage.getRegionIndex(xIndex,yIndex,0) == regionIndex){
                                numOctantsToRemove += (8 - totalOctants);
                            }
                            // lg.info("+z side including nothing");
                            if (regionInfo.isIndexInRegion(xIndex + yIndex * numX + (numZ - 1) * numXY)) {
                                // (regionImage.getRegionIndex(xIndex,yIndex,numZ-1) == regionIndex){
                                numOctantsToRemove += (8 - totalOctants);
                            }
                        }
                    }
                    size -= sizeOfPixel * 0.125 * numOctantsToRemove;
                    if (lg.isInfoEnabled()) {
                        lg.info("size=" + size);
                    }
                    break;
                }
        }
        VolumeGeometricRegion volumeRegion = new VolumeGeometricRegion(name, size, volumeUnit, subVolume, regionInfo.getRegionIndex());
        regionList.add(volumeRegion);
        if (lg.isInfoEnabled()) {
            lg.info("added volumeRegion(" + volumeRegion.getName() + ")");
        }
    }
    // 
    for (int i = 0; i < surfaceCollection.getSurfaceCount(); i++) {
        cbit.vcell.geometry.surface.Surface surface = surfaceCollection.getSurfaces(i);
        int exteriorRegionIndex = surface.getExteriorRegionIndex();
        int interiorRegionIndex = surface.getInteriorRegionIndex();
        cbit.vcell.geometry.SubVolume surfaceExteriorSubvolume = geometrySpec.getSubVolume(regionInfos[exteriorRegionIndex].getPixelValue());
        cbit.vcell.geometry.SubVolume surfaceInteriorSubvolume = geometrySpec.getSubVolume(regionInfos[interiorRegionIndex].getPixelValue());
        String name = "membrane_" + surfaceInteriorSubvolume.getName() + interiorRegionIndex + "_" + surfaceExteriorSubvolume.getName() + exteriorRegionIndex;
        double correctedArea = surface.getArea();
        if (geometrySpec.getDimension() == 2) {
            correctedArea /= geometrySpec.getExtent().getZ();
        } else if (geometrySpec.getDimension() == 1) {
            correctedArea /= (geometrySpec.getExtent().getY() * geometrySpec.getExtent().getZ());
        }
        SurfaceGeometricRegion surfaceRegion = new SurfaceGeometricRegion(name, correctedArea, surfaceUnit);
        regionList.add(surfaceRegion);
        // 
        // connect this surfaceLocation to its exterior and interior volumeLocations
        // 
        VolumeGeometricRegion exteriorVolumeRegion = null;
        for (int j = 0; j < regionList.size(); j++) {
            GeometricRegion region = regionList.elementAt(j);
            if (region instanceof VolumeGeometricRegion && region.getName().equals(surfaceExteriorSubvolume.getName() + exteriorRegionIndex)) {
                exteriorVolumeRegion = (VolumeGeometricRegion) region;
            }
        }
        surfaceRegion.addAdjacentGeometricRegion(exteriorVolumeRegion);
        exteriorVolumeRegion.addAdjacentGeometricRegion(surfaceRegion);
        VolumeGeometricRegion interiorVolumeRegion = null;
        for (int j = 0; j < regionList.size(); j++) {
            GeometricRegion region = regionList.elementAt(j);
            if (region instanceof VolumeGeometricRegion && region.getName().equals(surfaceInteriorSubvolume.getName() + interiorRegionIndex)) {
                interiorVolumeRegion = (VolumeGeometricRegion) region;
            }
        }
        surfaceRegion.addAdjacentGeometricRegion(interiorVolumeRegion);
        interiorVolumeRegion.addAdjacentGeometricRegion(surfaceRegion);
        if (lg.isInfoEnabled()) {
            lg.info("added surfaceRegion(" + surfaceRegion.getName() + ")");
        }
    }
    return org.vcell.util.BeanUtils.getArray(regionList, GeometricRegion.class);
}

Example 23 with RegionImage

use of cbit.vcell.geometry.RegionImage in project vcell by virtualcell.

the class GeometrySurfaceUtils method getUpdatedRegionImage.

/**
 * Insert the method's description here.
 * Creation date: (6/28/2004 2:56:03 PM)
 * @return cbit.vcell.geometry.RegionImage
 * @param geoSurfaceDescription cbit.vcell.geometry.surface.GeometrySurfaceDescription
 */
static cbit.vcell.geometry.RegionImage getUpdatedRegionImage(GeometrySurfaceDescription geoSurfaceDescription) throws cbit.vcell.geometry.GeometryException, cbit.image.ImageException, cbit.vcell.parser.ExpressionException {
    // 
    // make new RegionImage if necessary missing or wrong size
    // 
    cbit.vcell.geometry.GeometrySpec geometrySpec = geoSurfaceDescription.getGeometry().getGeometrySpec();
    if (geoSurfaceDescription.getRegionImage() == null || geoSurfaceDescription.getFilterCutoffFrequency() != geoSurfaceDescription.getRegionImage().getFilterCutoffFrequency() || geoSurfaceDescription.getRegionImage().getNumX() != geoSurfaceDescription.getVolumeSampleSize().getX() || geoSurfaceDescription.getRegionImage().getNumY() != geoSurfaceDescription.getVolumeSampleSize().getY() || geoSurfaceDescription.getRegionImage().getNumZ() != geoSurfaceDescription.getVolumeSampleSize().getZ()) {
        cbit.image.VCImage image = geometrySpec.createSampledImage(geoSurfaceDescription.getVolumeSampleSize());
        cbit.vcell.geometry.RegionImage regionImage = new RegionImage(image, geoSurfaceDescription.getGeometry().getDimension(), geoSurfaceDescription.getGeometry().getExtent(), geoSurfaceDescription.getGeometry().getOrigin(), geoSurfaceDescription.getFilterCutoffFrequency().doubleValue());
        return regionImage;
    }
    return geoSurfaceDescription.getRegionImage();
}

Example 24 with RegionImage

use of cbit.vcell.geometry.RegionImage in project vcell by virtualcell.

the class GeometrySurfaceUtils method updateGeometricRegions.

/**
 * Insert the method's description here.
 * Creation date: (5/26/2004 10:19:59 AM)
 */
public static void updateGeometricRegions(GeometrySurfaceDescription geoSurfaceDescription) throws cbit.vcell.geometry.GeometryException, cbit.image.ImageException, cbit.vcell.parser.ExpressionException, java.beans.PropertyVetoException {
    // 
    // updates if necessary: regionImage, surfaceCollection and resolvedLocations[]
    // 
    // assumes that volumeSampleSize and filterCutoffFrequency are already specified
    // 
    // 
    // 
    // make new RegionImage if necessary missing or wrong size
    // 
    boolean bChanged = false;
    cbit.vcell.geometry.RegionImage updatedRegionImage = getUpdatedRegionImage(geoSurfaceDescription);
    if (updatedRegionImage != geoSurfaceDescription.getRegionImage()) {
        // getUpdatedRegionImage returns same image if not changed
        bChanged = true;
    }
    // 
    // parse regionImage into VolumeGeometricRegions and SurfaceCollection into SurfaceGeometricRegions
    // 
    GeometricRegion[] updatedGeometricRegions = geoSurfaceDescription.getGeometricRegions();
    if (geoSurfaceDescription.getGeometricRegions() == null || bChanged) {
        updatedGeometricRegions = getUpdatedGeometricRegions(geoSurfaceDescription, updatedRegionImage, updatedRegionImage.getSurfacecollection());
    }
    geoSurfaceDescription.setGeometricRegions(updatedGeometricRegions);
}

Example 25 with RegionImage

use of cbit.vcell.geometry.RegionImage in project vcell by virtualcell.

the class GeometryFileWriter method write.

/**
 * Insert the method's description here.
 * Creation date: (7/19/2004 10:54:30 AM)
 * @param geometrySurfaceDescription cbit.vcell.geometry.surface.GeometrySurfaceDescription
 * @throws IOException
 */
public static void write(Writer writer, Geometry resampledGeometry) throws IOException {
    // 
    // "name" name
    // "dimension" dimension
    // "extent" extentx extenty extentz
    // "origin" originx originy originz
    // "volumeRegions" num
    // name totalVolume featureHandle
    // "membraneRegions" num
    // name totalArea volumeRegionIndex1 volumeRegionIndex2
    // "volumeSamples" numX, numY, numZ
    // uncompressed regionIndexs for each volume element
    // compressed regionIndexs for each volume element
    // "nodes" num
    // nodeIndex x y z
    // "cells" num
    // cellIndex patchIndex node1 node2 node3 node4
    // "celldata"
    // insideVolumeIndex outsideVolumeIndex area normalx normaly normalz
    // 
    // 
    // When we are writing volume regions, we sort regions so that ID is equal to index
    // 
    writer.write("name " + resampledGeometry.getName() + "\n");
    writer.write("dimension " + resampledGeometry.getDimension() + "\n");
    org.vcell.util.Extent extent = resampledGeometry.getExtent();
    org.vcell.util.Origin origin = resampledGeometry.getOrigin();
    switch(resampledGeometry.getDimension()) {
        case 1:
            writer.write("size " + extent.getX() + "\n");
            writer.write("origin " + origin.getX() + "\n");
            break;
        case 2:
            writer.write("size " + extent.getX() + " " + extent.getY() + "\n");
            writer.write("origin " + origin.getX() + " " + origin.getY() + "\n");
            break;
        case 3:
            writer.write("size " + extent.getX() + " " + extent.getY() + " " + extent.getZ() + "\n");
            writer.write("origin " + origin.getX() + " " + origin.getY() + " " + origin.getZ() + "\n");
            break;
    }
    GeometrySurfaceDescription geoSurfaceDesc = resampledGeometry.getGeometrySurfaceDescription();
    RegionImage regionImage = geoSurfaceDesc.getRegionImage();
    SurfaceCollection surfaceCollection = geoSurfaceDesc.getSurfaceCollection();
    GeometricRegion[] geometricRegions = geoSurfaceDesc.getGeometricRegions();
    int numVolumeRegions = 0;
    int numMembraneRegions = 0;
    Vector<VolumeGeometricRegion> volRegionList = new Vector<VolumeGeometricRegion>();
    if (geometricRegions != null) {
        for (int i = 0; i < geometricRegions.length; i++) {
            if (geometricRegions[i] instanceof VolumeGeometricRegion) {
                numVolumeRegions++;
                volRegionList.add((VolumeGeometricRegion) geometricRegions[i]);
            } else if (geometricRegions[i] instanceof SurfaceGeometricRegion) {
                numMembraneRegions++;
            }
        }
    }
    // 
    // get ordered array of volume regions (where "id" == index into array)... fail if impossible
    // 
    java.util.Collections.sort(volRegionList, new Comparator<VolumeGeometricRegion>() {

        public int compare(VolumeGeometricRegion reg1, VolumeGeometricRegion reg2) {
            if (reg1.getRegionID() < reg2.getRegionID()) {
                return -1;
            } else if (reg1.getRegionID() > reg2.getRegionID()) {
                return 1;
            } else {
                return 0;
            }
        }

        public boolean equals(Object obj) {
            return this == obj;
        }
    });
    VolumeGeometricRegion[] volRegions = (VolumeGeometricRegion[]) org.vcell.util.BeanUtils.getArray(volRegionList, VolumeGeometricRegion.class);
    writer.write("volumeRegions " + numVolumeRegions + "\n");
    for (int i = 0; i < volRegions.length; i++) {
        if (volRegions[i].getRegionID() != i) {
            throw new RuntimeException("Region ID != Region Index, they must be the same!");
        }
        writer.write(volRegions[i].getName() + " " + volRegions[i].getSize() + " " + volRegions[i].getSubVolume().getHandle() + "\n");
    }
    writer.write("membraneRegions " + numMembraneRegions + "\n");
    if (geometricRegions != null) {
        for (int i = 0; i < geometricRegions.length; i++) {
            if (geometricRegions[i] instanceof SurfaceGeometricRegion) {
                SurfaceGeometricRegion surfaceRegion = (SurfaceGeometricRegion) geometricRegions[i];
                GeometricRegion[] neighbors = surfaceRegion.getAdjacentGeometricRegions();
                VolumeGeometricRegion insideRegion = (VolumeGeometricRegion) neighbors[0];
                VolumeGeometricRegion outsideRegion = (VolumeGeometricRegion) neighbors[1];
                writer.write(surfaceRegion.getName() + " " + surfaceRegion.getSize() + " " + insideRegion.getRegionID() + " " + outsideRegion.getRegionID() + "\n");
            }
        }
    }
    // 
    // write volume samples
    // 
    ISize volumeSampleSize = geoSurfaceDesc.getVolumeSampleSize();
    switch(resampledGeometry.getDimension()) {
        case 1:
            writer.write("volumeSamples " + volumeSampleSize.getX() + "\n");
            break;
        case 2:
            writer.write("volumeSamples " + volumeSampleSize.getX() + " " + volumeSampleSize.getY() + "\n");
            break;
        case 3:
            writer.write("volumeSamples " + volumeSampleSize.getX() + " " + volumeSampleSize.getY() + " " + volumeSampleSize.getZ() + "\n");
            break;
    }
    // regionImage
    if (regionImage != null) {
        if (regionImage.getNumRegions() > 65536) {
            throw new RuntimeException("cannot process a geometry with more than 65536 volume regions");
        }
        byte[] uncompressedRegionIDs = new byte[2 * regionImage.getNumX() * regionImage.getNumY() * regionImage.getNumZ()];
        for (int i = 0, j = 0; i < uncompressedRegionIDs.length; i += 2, j++) {
            int regindex = regionImage.getRegionInfoFromOffset(j).getRegionIndex();
            uncompressedRegionIDs[i] = (byte) (regindex & 0x000000ff);
            uncompressedRegionIDs[i + 1] = (byte) ((regindex & 0x0000ff00) >> 8);
        }
        ByteArrayOutputStream bos = new ByteArrayOutputStream();
        DeflaterOutputStream dos = new DeflaterOutputStream(bos);
        dos.write(uncompressedRegionIDs, 0, uncompressedRegionIDs.length);
        dos.close();
        byte[] compressedRegionIDs = bos.toByteArray();
        writer.write(org.vcell.util.Hex.toString(compressedRegionIDs) + "\n");
    } else {
        writer.write("\n");
    }
    // 
    if (surfaceCollection == null) {
        throw new RuntimeException("geometry is not updated");
    }
    int numCells = surfaceCollection.getTotalPolygonCount();
    writer.write("cells " + numCells + "\n");
    // "celldata"
    // insideVolumeIndex outsideVolumeIndex area normalx normaly normalz
    // 
    int cellID = 0;
    int dimension = resampledGeometry.getDimension();
    double correctCoeff = 1;
    if (dimension == 1) {
        correctCoeff = extent.getY() * extent.getZ();
    } else if (dimension == 2) {
        correctCoeff = extent.getZ();
    }
    if (surfaceCollection != null) {
        for (int i = 0; i < surfaceCollection.getSurfaceCount(); i++) {
            Surface surface = surfaceCollection.getSurfaces(i);
            int region1Outside = 0;
            int region1Inside = 0;
            for (int j = 0; j < surface.getPolygonCount(); j++) {
                Quadrilateral polygon = (Quadrilateral) surface.getPolygons(j);
                Node[] node = polygon.getNodes();
                cbit.vcell.render.Vect3d elementCoord = new cbit.vcell.render.Vect3d();
                int nodesOnBoundary = 0;
                for (int k = 0; k < node.length; k++) {
                    if (!node[k].getMoveX() || (dimension > 1 && !node[k].getMoveY()) || (dimension == 3 && !node[k].getMoveZ())) {
                        nodesOnBoundary++;
                    }
                }
                if (nodesOnBoundary == 0) {
                    for (int k = 0; k < node.length; k++) {
                        elementCoord.add(new cbit.vcell.render.Vect3d(node[k].getX(), node[k].getY(), node[k].getZ()));
                    }
                    elementCoord.scale(0.25);
                } else if (nodesOnBoundary == 2) {
                    for (int k = 0; k < node.length; k++) {
                        if (!node[k].getMoveX() || !node[k].getMoveY() || !node[k].getMoveZ()) {
                            elementCoord.add(new cbit.vcell.render.Vect3d(node[k].getX(), node[k].getY(), node[k].getZ()));
                        }
                    }
                    elementCoord.scale(0.5);
                } else if (nodesOnBoundary == 3) {
                    for (int k = 0; k < node.length; k++) {
                        if (!node[k].getMoveX() && !node[k].getMoveY() || !node[k].getMoveY() && !node[k].getMoveZ() || !node[k].getMoveX() && !node[k].getMoveZ()) {
                            elementCoord.set(node[k].getX(), node[k].getY(), node[k].getZ());
                        }
                    }
                } else {
                    throw new RuntimeException("Unexcepted number of nodes on boundary for a polygon: " + nodesOnBoundary);
                }
                cbit.vcell.render.Vect3d unitNormal = new cbit.vcell.render.Vect3d();
                polygon.getUnitNormal(unitNormal);
                int volNeighbor1Region = regionImage.getRegionInfoFromOffset(polygon.getVolIndexNeighbor1()).getRegionIndex();
                int volNeighbor2Region = regionImage.getRegionInfoFromOffset(polygon.getVolIndexNeighbor2()).getRegionIndex();
                if (surface.getExteriorRegionIndex() == volNeighbor1Region && surface.getInteriorRegionIndex() == volNeighbor2Region) {
                    region1Outside++;
                }
                if (surface.getExteriorRegionIndex() == volNeighbor2Region && surface.getInteriorRegionIndex() == volNeighbor1Region) {
                    region1Inside++;
                }
                writer.write(cellID + " " + polygon.getVolIndexNeighbor1() + " " + polygon.getVolIndexNeighbor2() + " " + polygon.getArea() / correctCoeff + " " + elementCoord.getX() + " " + elementCoord.getY() + " " + elementCoord.getZ() + " " + unitNormal.getX() + " " + unitNormal.getY() + " " + unitNormal.getZ() + "\n");
                cellID++;
            }
            if (region1Inside != surface.getPolygonCount() && region1Outside != surface.getPolygonCount()) {
                throw new RuntimeException("Volume neighbor regions not consistent: [total, inside, outside]=" + surface.getPolygonCount() + "," + region1Inside + "," + region1Outside + "]");
            }
        }
    }
}