Examples with Origin org.vcell.util.Origin used on opensource projects

Example 11 with Origin

use of org.vcell.util.Origin in project vcell by virtualcell.

the class CartesianMesh method read.

/**
 * This method was created by a SmartGuide.
 * @param tokens java.util.StringTokenizer
 * @exception java.lang.Exception The exception description.
 */
private void read(CommentStringTokenizer tokens, MembraneMeshMetrics membraneMeshMetrics) throws MathException {
    // 
    // clear previous contents
    // 
    membraneElements = null;
    // 
    // read new stuff
    // 
    String token = null;
    token = tokens.nextToken();
    if (token.equalsIgnoreCase(VCML.Version)) {
        // 
        // read version number
        // 
        token = tokens.nextToken();
        this.version = token;
        token = tokens.nextToken();
    }
    if (token.equalsIgnoreCase(VCML.CartesianMesh)) {
        token = tokens.nextToken();
    } else {
        throw new MathFormatException("unexpected token " + token + " expecting " + VCML.CartesianMesh);
    }
    // 
    // only Version 1.1 and later supports membrane connectivity  (as of 8/30/2000)
    // 
    boolean bConnectivity = false;
    if (version.equals(VERSION_1_1) || version.equals(VERSION_1_2)) {
        bConnectivity = true;
    }
    // 
    // only Version 1.2 and later supports Regions
    // 
    boolean bRegions = false;
    if (version.equals(VERSION_1_2)) {
        bRegions = true;
        meshRegionInfo = new MeshRegionInfo();
    }
    if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
        throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
    }
    while (tokens.hasMoreTokens()) {
        token = tokens.nextToken();
        if (token.equalsIgnoreCase(VCML.EndBlock)) {
            break;
        }
        if (token.equalsIgnoreCase(VCML.Size)) {
            int sx, sy, sz;
            try {
                token = tokens.nextToken();
                sx = Integer.valueOf(token).intValue();
                token = tokens.nextToken();
                sy = Integer.valueOf(token).intValue();
                token = tokens.nextToken();
                sz = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("expected:  " + VCML.Size + " # # #");
            }
            setSize(sx, sy, sz);
            continue;
        }
        if (token.equalsIgnoreCase(VCML.Extent)) {
            double ex, ey, ez;
            try {
                token = tokens.nextToken();
                ex = Double.valueOf(token).doubleValue();
                token = tokens.nextToken();
                ey = Double.valueOf(token).doubleValue();
                token = tokens.nextToken();
                ez = Double.valueOf(token).doubleValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("expected:  " + VCML.Extent + " # # #");
            }
            setExtent(new Extent(ex, ey, ez));
            continue;
        }
        if (token.equalsIgnoreCase(VCML.Origin)) {
            double ox, oy, oz;
            try {
                token = tokens.nextToken();
                ox = Double.valueOf(token).doubleValue();
                token = tokens.nextToken();
                oy = Double.valueOf(token).doubleValue();
                token = tokens.nextToken();
                oz = Double.valueOf(token).doubleValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("expected:  " + VCML.Origin + " # # #");
            }
            setOrigin(new Origin(ox, oy, oz));
            continue;
        }
        // 
        if (token.equalsIgnoreCase(VCML.VolumeRegionsMapSubvolume)) {
            token = tokens.nextToken();
            if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
                throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
            }
            token = tokens.nextToken();
            int numVolumeRegions = 0;
            try {
                numVolumeRegions = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("unexpected token " + token + " expecting the VolumeRegionsMapSubvolume list length");
            }
            int checkCount = 0;
            while (tokens.hasMoreTokens()) {
                token = tokens.nextToken();
                if (token.equalsIgnoreCase(VCML.EndBlock)) {
                    break;
                }
                try {
                    int volRegionID = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    int subvolumeID = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    double volume = Double.valueOf(token).doubleValue();
                    String subdomainName = null;
                    if (subdomainInfo != null) {
                        subdomainName = subdomainInfo.getCompartmentSubdomainName(subvolumeID);
                    }
                    meshRegionInfo.mapVolumeRegionToSubvolume(volRegionID, subvolumeID, volume, subdomainName);
                } catch (NumberFormatException e) {
                    throw new MathFormatException("expected:  # # #");
                }
                checkCount += 1;
            }
            if (checkCount != numVolumeRegions) {
                throw new MathFormatException("CartesianMesh.read->VolumeRegionsMapSubvolume: read " + checkCount + " VolRegions but was expecting " + numVolumeRegions);
            }
            continue;
        }
        if (token.equalsIgnoreCase(VCML.MembraneRegionsMapVolumeRegion)) {
            token = tokens.nextToken();
            if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
                throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
            }
            token = tokens.nextToken();
            int numMembraneRegions = 0;
            try {
                numMembraneRegions = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("unexpected token " + token + " expecting the MembraneRegionsMapVolumeRegion list length");
            }
            int checkCount = 0;
            while (tokens.hasMoreTokens()) {
                token = tokens.nextToken();
                if (token.equalsIgnoreCase(VCML.EndBlock)) {
                    break;
                }
                try {
                    int memRegionID = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    int volRegionIn = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    int volRegionOut = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    double surface = Double.valueOf(token).doubleValue();
                    meshRegionInfo.mapMembraneRegionToVolumeRegion(memRegionID, volRegionIn, volRegionOut, surface);
                } catch (NumberFormatException e) {
                    throw new MathFormatException("expected:  # # #");
                }
                checkCount += 1;
            }
            if (checkCount != numMembraneRegions) {
                throw new MathFormatException("CartesianMesh.read->MembraneRegionsMapVolumeRegion: read " + checkCount + " MembraneRegions but was expecting " + numMembraneRegions);
            }
            continue;
        }
        if (token.equalsIgnoreCase(VCML.VolumeElementsMapVolumeRegion)) {
            token = tokens.nextToken();
            if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
                throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
            }
            token = tokens.nextToken();
            int numVolumeElements = 0;
            try {
                numVolumeElements = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("unexpected token " + token + " expecting the VolumeElementsMapVolumeRegion list length");
            }
            token = tokens.nextToken();
            boolean bCompressed = token.equalsIgnoreCase("Compressed");
            if (!bCompressed) {
                if (!token.equalsIgnoreCase("UnCompressed")) {
                    throw new MathFormatException("unexpected token " + token + " expecting Compress or UnCompress");
                }
            }
            byte[] volumeElementMap = new byte[numVolumeElements];
            int checkCount = 0;
            if (bCompressed) {
                // Get HEX encoded bytes of the compressed VolumeElements-RegionID Map
                StringBuffer hexOfCompressed = new StringBuffer();
                while (tokens.hasMoreTokens()) {
                    token = tokens.nextToken();
                    if (token.equalsIgnoreCase(VCML.EndBlock)) {
                        break;
                    }
                    hexOfCompressed.append(token);
                }
                // Un-HEX the compressed data
                byte[] compressedData = Hex.toBytes(hexOfCompressed.toString());
                try {
                    meshRegionInfo.setCompressedVolumeElementMapVolumeRegion(compressedData, numVolumeElements);
                } catch (IOException e) {
                    throw new MathFormatException("CartesianMesh.read->VolumeElementsMapVolumeRegion " + e.toString());
                }
                checkCount = meshRegionInfo.getUncompressedVolumeElementMapVolumeRegionLength();
            } else {
                while (tokens.hasMoreTokens()) {
                    token = tokens.nextToken();
                    if (token.equalsIgnoreCase(VCML.EndBlock)) {
                        break;
                    }
                    try {
                        int volumeRegionID = Integer.valueOf(token).intValue();
                        volumeElementMap[checkCount] = (byte) volumeRegionID;
                    } catch (NumberFormatException e) {
                        throw new MathFormatException("expected:  # # #");
                    }
                    checkCount += 1;
                }
            }
            if (checkCount != numVolumeElements && checkCount != 2 * numVolumeElements) {
                throw new MathFormatException("CartesianMesh.read->VolumeElementsMapVolumeRegion: read " + checkCount + " VolumeElements but was expecting " + numVolumeElements);
            }
            continue;
        }
        // 
        // 
        // 
        HashMap<Integer, Integer> volumeRegionMapSubvolume = getVolumeRegionMapSubvolume();
        if (token.equalsIgnoreCase(VCML.MembraneElements)) {
            // 
            // read '{'
            // 
            token = tokens.nextToken();
            if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
                throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
            }
            token = tokens.nextToken();
            int numMemElements = 0;
            try {
                numMemElements = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("unexpected token " + token + " expecting the membraneElement list length");
            }
            // 
            // read list of the following format:
            // 
            // memIndex insideVolIndex outsideVolIndex
            // 
            membraneElements = new MembraneElement[numMemElements];
            int index = 0;
            int[] membraneElementMapMembraneRegion = null;
            if (bRegions) {
                membraneElementMapMembraneRegion = new int[numMemElements];
                meshRegionInfo.mapMembraneElementsToMembraneRegions(membraneElementMapMembraneRegion);
            }
            // 
            while (tokens.hasMoreTokens()) {
                token = tokens.nextToken();
                if (token.equalsIgnoreCase(VCML.EndBlock)) {
                    break;
                }
                int memIndex = -1;
                int insideIndex = -1;
                int outsideIndex = -1;
                try {
                    // 
                    // read first three tokens of a membrane element
                    // 
                    // membraneIndex   insideIndex    outsideIndex
                    // 
                    memIndex = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    insideIndex = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    outsideIndex = Integer.valueOf(token).intValue();
                    if (subdomainInfo != null) {
                        int insideRegionIndex = getVolumeRegionIndex(insideIndex);
                        int outsideRegionIndex = getVolumeRegionIndex(outsideIndex);
                        int insideSubVolumeHandle = volumeRegionMapSubvolume.get(insideRegionIndex);
                        int outsideSubVolumeHandle = volumeRegionMapSubvolume.get(outsideRegionIndex);
                        int realInsideSubVolumeHandle = subdomainInfo.getInside(insideSubVolumeHandle, outsideSubVolumeHandle);
                        if (realInsideSubVolumeHandle != insideSubVolumeHandle) {
                            int temp = insideIndex;
                            insideIndex = outsideIndex;
                            outsideIndex = temp;
                        }
                    }
                } catch (NumberFormatException e) {
                    throw new MathFormatException("expected:  # # #");
                }
                MembraneElement me = null;
                // 
                if (bConnectivity) {
                    try {
                        token = tokens.nextToken();
                        int neighbor1 = Integer.valueOf(token).intValue();
                        token = tokens.nextToken();
                        int neighbor2 = Integer.valueOf(token).intValue();
                        token = tokens.nextToken();
                        int neighbor3 = Integer.valueOf(token).intValue();
                        token = tokens.nextToken();
                        int neighbor4 = Integer.valueOf(token).intValue();
                        // 
                        if (bRegions) {
                            token = tokens.nextToken();
                            int regionID = Integer.valueOf(token).intValue();
                            membraneElementMapMembraneRegion[memIndex] = regionID;
                        }
                        if (membraneMeshMetrics == null) {
                            me = new MembraneElement(memIndex, insideIndex, outsideIndex, neighbor1, neighbor2, neighbor3, neighbor4, MembraneElement.AREA_UNDEFINED, 0, 0, 0, 0, 0, 0);
                        } else {
                            me = new MembraneElement(memIndex, insideIndex, outsideIndex, neighbor1, neighbor2, neighbor3, neighbor4, membraneMeshMetrics.areas[memIndex], membraneMeshMetrics.normals[memIndex][0], membraneMeshMetrics.normals[memIndex][1], membraneMeshMetrics.normals[memIndex][2], membraneMeshMetrics.centroids[memIndex][0], membraneMeshMetrics.centroids[memIndex][1], membraneMeshMetrics.centroids[memIndex][2]);
                        }
                    } catch (NumberFormatException e) {
                        throw new MathFormatException("expected:  # # # # # # #");
                    }
                } else {
                    me = new MembraneElement(memIndex, insideIndex, outsideIndex);
                }
                membraneElements[index] = me;
                index++;
            }
            continue;
        }
        if (token.equalsIgnoreCase(VCML.ContourElements)) {
            // 
            // read '{'
            // 
            token = tokens.nextToken();
            if (!token.equalsIgnoreCase(VCML.BeginBlock)) {
                throw new MathFormatException("unexpected token " + token + " expecting " + VCML.BeginBlock);
            }
            token = tokens.nextToken();
            int numContourElements = 0;
            try {
                numContourElements = Integer.valueOf(token).intValue();
            } catch (NumberFormatException e) {
                throw new MathFormatException("unexpected token " + token + " expecting the contourElement list length");
            }
            // 
            // read list of the following format:
            // 
            // contourIndex volumeIndex beginCoord endCoord prevIndex nextIndex
            // 
            contourElements = new ContourElement[numContourElements];
            int index = 0;
            // 
            while (tokens.hasMoreTokens()) {
                token = tokens.nextToken();
                if (token.equalsIgnoreCase(VCML.EndBlock)) {
                    break;
                }
                ContourElement ce = null;
                try {
                    // 
                    // read first two tokens of a contour element
                    // 
                    // contourIndex volumeIndex
                    // 
                    int contourIndex = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    int volumeIndex = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    // 
                    // read beginCoord endCoord
                    // 
                    double beginX = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    double beginY = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    double beginZ = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    double endX = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    double endY = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    double endZ = Double.valueOf(token).doubleValue();
                    token = tokens.nextToken();
                    Coordinate begin = new Coordinate(beginX, beginY, beginZ);
                    Coordinate end = new Coordinate(endX, endY, endZ);
                    // 
                    // read last two tokens of a contour element
                    // 
                    // prevContourIndex nextContourIndex
                    // 
                    int prevContourIndex = Integer.valueOf(token).intValue();
                    token = tokens.nextToken();
                    int nextContourIndex = Integer.valueOf(token).intValue();
                    ce = new ContourElement(contourIndex, volumeIndex, begin, end, prevContourIndex, nextContourIndex);
                } catch (NumberFormatException e) {
                    throw new MathFormatException("expected:  %d %d   %f %f %f    %f %f %f   %d %d");
                }
                contourElements[index] = ce;
                index++;
            }
            continue;
        }
        throw new MathFormatException("unexpected identifier " + token);
    }
    switch(getGeometryDimension()) {
        case 1:
            {
                if (extent.getY() != 1 || extent.getZ() != 1) {
                    System.out.println("Extent " + extent.toString() + " for a 1-D mesh truncated to 1 for y and z");
                    setExtent(new Extent(extent.getX(), 1.0, 1.0));
                }
                break;
            }
        case 2:
            {
                if (extent.getZ() != 1) {
                    System.out.println("Extent " + extent.toString() + " for a 2-D mesh truncated to 1 for z");
                    setExtent(new Extent(extent.getX(), extent.getY(), 1.0));
                }
                break;
            }
    }
}

Example 12 with Origin

use of org.vcell.util.Origin in project vcell by virtualcell.

the class CartesianMeshMovingBoundary method readMeshFile.

public static CartesianMeshMovingBoundary readMeshFile(File meshFile) throws Exception {
    CartesianMeshMovingBoundary mesh = new CartesianMeshMovingBoundary();
    if (H5.H5open() < 0) {
        throw new Exception("H5.H5open() failed");
    }
    FileFormat fileFormat = FileFormat.getFileFormat(FileFormat.FILE_TYPE_HDF5);
    if (fileFormat == null) {
        throw new Exception("FileFormat.getFileFormat(FileFormat.FILE_TYPE_HDF5) failed, returned null.");
    }
    FileFormat meshH5File = null;
    try {
        meshH5File = fileFormat.createInstance(meshFile.getAbsolutePath(), FileFormat.READ);
        meshH5File.open();
        DefaultMutableTreeNode rootNode = (DefaultMutableTreeNode) meshH5File.getRootNode();
        Group rootGroup = (Group) rootNode.getUserObject();
        Group meshGroup = null;
        for (Object member : rootGroup.getMemberList()) {
            if (member instanceof Group) {
                Group g = (Group) member;
                if (g.getName().equals(Group_Mesh))
                    ;
                {
                    meshGroup = g;
                    break;
                }
            }
        }
        if (meshGroup == null) {
            throw new Exception(Group_Mesh + " group not found in mesh");
        }
        for (Object member : meshGroup.getMemberList()) {
            if (member instanceof Dataset) {
                Dataset ds = (Dataset) member;
                Object data = ds.getData();
                MeshDataset mds = MeshDataset.valueOf(ds.getName());
                switch(mds) {
                    case dimension:
                        mesh.dimension = ((int[]) data)[0];
                        break;
                    case extent:
                        {
                            double[] darr = (double[]) data;
                            mesh.extent = new Extent(darr[0], darr[1], 0.5);
                            break;
                        }
                    case origin:
                        {
                            double[] darr = (double[]) data;
                            mesh.origin = new Origin(darr[0], darr[1], 0.5);
                            break;
                        }
                    case size:
                        {
                            int[] iarr = (int[]) data;
                            mesh.size = new ISize(iarr[0], iarr[1], 1);
                            break;
                        }
                }
            }
        }
    } finally {
        if (meshH5File != null) {
            meshH5File.close();
        }
    }
    return mesh;
}

Example 13 with Origin

use of org.vcell.util.Origin in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeChomboSpec.

private void writeChomboSpec() throws ExpressionException, SolverException, PropertyVetoException, ClassNotFoundException, IOException, GeometryException, ImageException {
    if (!bChomboSolver) {
        return;
    }
    GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
    int dimension = geometrySpec.getDimension();
    if (dimension == 1) {
        throw new SolverException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " is only supported for simulations with 2D or 3D geometry.");
    }
    Simulation simulation = getSimulationTask().getSimulation();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    ChomboSolverSpec chomboSolverSpec = solverTaskDescription.getChomboSolverSpec();
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_BEGIN);
    printWriter.println(FVInputFileKeyword.DIMENSION + " " + geometrySpec.getDimension());
    Extent extent = geometrySpec.getExtent();
    Origin origin = geometrySpec.getOrigin();
    ISize isize = simulation.getMeshSpecification().getSamplingSize();
    switch(geometrySpec.getDimension()) {
        case 2:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY());
            break;
        case 3:
            printWriter.println(FVInputFileKeyword.MESH_SIZE + " " + isize.getX() + " " + isize.getY() + " " + isize.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_SIZE + " " + extent.getX() + " " + extent.getY() + " " + extent.getZ());
            printWriter.println(FVInputFileKeyword.DOMAIN_ORIGIN + " " + origin.getX() + " " + origin.getY() + " " + origin.getZ());
            break;
    }
    List<CompartmentSubDomain> featureList = new ArrayList<CompartmentSubDomain>();
    Enumeration<SubDomain> enum1 = simulation.getMathDescription().getSubDomains();
    while (enum1.hasMoreElements()) {
        SubDomain sd = enum1.nextElement();
        if (sd instanceof CompartmentSubDomain) {
            featureList.add((CompartmentSubDomain) sd);
        }
    }
    int numFeatures = featureList.size();
    CompartmentSubDomain[] features = featureList.toArray(new CompartmentSubDomain[0]);
    int[] phases = new int[numFeatures];
    Arrays.fill(phases, -1);
    phases[numFeatures - 1] = 0;
    int[] numAssigned = new int[] { 1 };
    assignPhases(features, numFeatures - 1, phases, numAssigned);
    Map<String, Integer> subDomainPhaseMap = new HashMap<String, Integer>();
    for (int i = 0; i < phases.length; ++i) {
        if (phases[i] == -1) {
            throw new SolverException("Failed to assign a phase to CompartmentSubdomain '" + features[i].getName() + "'. It might be caused by too coarsh a mesh.");
        }
        subDomainPhaseMap.put(features[i].getName(), phases[i]);
    }
    SubVolume[] subVolumes = geometrySpec.getSubVolumes();
    if (geometrySpec.hasImage()) {
        Geometry geometry = (Geometry) BeanUtils.cloneSerializable(simulation.getMathDescription().getGeometry());
        Geometry simGeometry = geometry;
        VCImage img = geometry.getGeometrySpec().getImage();
        int factor = Math.max(Math.max(img.getNumX(), img.getNumY()), img.getNumZ()) < 512 ? 2 : 1;
        ISize distanceMapMeshSize = new ISize(img.getNumX() * factor, img.getNumY() * factor, img.getNumZ() * factor);
        Vect3d deltaX = null;
        boolean bCellCentered = false;
        double dx = 0.5;
        double dy = 0.5;
        double dz = 0.5;
        int Nx = distanceMapMeshSize.getX();
        int Ny = distanceMapMeshSize.getY();
        int Nz = distanceMapMeshSize.getZ();
        if (dimension == 2) {
            // pad the 2D image with itself in order to obtain a 3D image used to compute the distance map
            // because the distance map algorithm is 3D only (using distance to triangles)
            byte[] oldPixels = img.getPixels();
            byte[] newPixels = new byte[oldPixels.length * 3];
            System.arraycopy(oldPixels, 0, newPixels, 0, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length, oldPixels.length);
            System.arraycopy(oldPixels, 0, newPixels, oldPixels.length * 2, oldPixels.length);
            double distX = geometry.getExtent().getX() / img.getNumX();
            double distY = geometry.getExtent().getY() / img.getNumY();
            // we set the distance on the z axis to something that makes sense
            double distZ = Math.max(distX, distY);
            Extent newExtent = new Extent(geometry.getExtent().getX(), geometry.getExtent().getY(), distZ * 3);
            VCImage newImage = new VCImageUncompressed(null, newPixels, newExtent, img.getNumX(), img.getNumY(), 3);
            // copy the pixel classes too
            ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
            for (VCPixelClass origPixelClass : geometry.getGeometrySpec().getImage().getPixelClasses()) {
                SubVolume origSubvolume = geometry.getGeometrySpec().getImageSubVolumeFromPixelValue(origPixelClass.getPixel());
                newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origPixelClass.getPixel()));
            }
            newImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
            simGeometry = new Geometry(geometry, newImage);
            Nz = 3;
        }
        GeometrySpec simGeometrySpec = simGeometry.getGeometrySpec();
        Extent simExtent = simGeometrySpec.getExtent();
        dx = simExtent.getX() / (Nx - 1);
        dy = simExtent.getY() / (Ny - 1);
        dz = simExtent.getZ() / (Nz - 1);
        if (Math.abs(dx - dy) > 0.1 * Math.max(dx, dy)) {
            dx = Math.min(dx, dy);
            dy = dx;
            Nx = (int) (simExtent.getX() / dx + 1);
            Ny = (int) (simExtent.getY() / dx + 1);
            if (dimension == 3) {
                dz = dx;
                Nz = (int) (simExtent.getZ() / dx + 1);
            }
        }
        deltaX = new Vect3d(dx, dy, dz);
        // one more point in each direction
        distanceMapMeshSize = new ISize(Nx + 1, Ny + 1, Nz + 1);
        Extent distanceMapExtent = new Extent(simExtent.getX() + dx, simExtent.getY() + dy, simExtent.getZ() + dz);
        simGeometrySpec.setExtent(distanceMapExtent);
        GeometrySurfaceDescription geoSurfaceDesc = simGeometry.getGeometrySurfaceDescription();
        geoSurfaceDesc.setVolumeSampleSize(distanceMapMeshSize);
        geoSurfaceDesc.updateAll();
        VCImage vcImage = RayCaster.sampleGeometry(simGeometry, distanceMapMeshSize, bCellCentered);
        SubvolumeSignedDistanceMap[] distanceMaps = DistanceMapGenerator.computeDistanceMaps(simGeometry, vcImage, bCellCentered);
        if (dimension == 2) {
            distanceMaps = DistanceMapGenerator.extractMiddleSlice(distanceMaps);
        }
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + simGeometrySpec.getNumSubVolumes() + " " + FVInputFileKeyword.DISTANCE_MAP);
        for (int i = 0; i < subVolumes.length; i++) {
            File distanceMapFile = new File(workingDirectory, getSimulationTask().getSimulationJobID() + "_" + subVolumes[i].getName() + DISTANCE_MAP_FILE_EXTENSION);
            writeDistanceMapFile(deltaX, distanceMaps[i], distanceMapFile);
            int phase = subDomainPhaseMap.get(subVolumes[i].getName());
            printWriter.println(subVolumes[i].getName() + " " + phase + " " + distanceMapFile.getAbsolutePath());
        }
    } else {
        printWriter.println(FVInputFileKeyword.SUBDOMAINS + " " + geometrySpec.getNumSubVolumes());
        Expression[] rvachevExps = convertAnalyticGeometryToRvachevFunction(geometrySpec);
        for (int i = 0; i < subVolumes.length; i++) {
            if (subVolumes[i] instanceof AnalyticSubVolume) {
                String name = subVolumes[i].getName();
                int phase = subDomainPhaseMap.get(name);
                printWriter.println(name + " " + phase + " ");
                printWriter.println(FVInputFileKeyword.IF + " " + rvachevExps[i].infix() + ";");
                printWriter.println(FVInputFileKeyword.USER + " " + ((AnalyticSubVolume) subVolumes[i]).getExpression().infix() + ";");
            }
        }
    }
    printWriter.println(FVInputFileKeyword.MAX_BOX_SIZE + " " + chomboSolverSpec.getMaxBoxSize());
    printWriter.println(FVInputFileKeyword.FILL_RATIO + " " + chomboSolverSpec.getFillRatio());
    printWriter.println(FVInputFileKeyword.RELATIVE_TOLERANCE + " " + simulation.getSolverTaskDescription().getErrorTolerance().getRelativeErrorTolerance());
    printWriter.println(FVInputFileKeyword.SAVE_VCELL_OUTPUT + " " + chomboSolverSpec.isSaveVCellOutput());
    printWriter.println(FVInputFileKeyword.SAVE_CHOMBO_OUTPUT + " " + chomboSolverSpec.isSaveChomboOutput());
    printWriter.println(FVInputFileKeyword.ACTIVATE_FEATURE_UNDER_DEVELOPMENT + " " + chomboSolverSpec.isActivateFeatureUnderDevelopment());
    printWriter.println(FVInputFileKeyword.SMALL_VOLFRAC_THRESHOLD + " " + chomboSolverSpec.getSmallVolfracThreshold());
    printWriter.println(FVInputFileKeyword.BLOCK_FACTOR + " " + chomboSolverSpec.getBlockFactor());
    printWriter.println(FVInputFileKeyword.TAGS_GROW + " " + chomboSolverSpec.getTagsGrow());
    // Refinement
    int numLevels = chomboSolverSpec.getNumRefinementLevels();
    // Refinements #Levels ratio 1, ratio 2, etc
    printWriter.print(FVInputFileKeyword.REFINEMENTS + " " + (numLevels + 1));
    List<Integer> ratios = chomboSolverSpec.getRefineRatioList();
    for (int i : ratios) {
        printWriter.print(" " + i);
    }
    // write last refinement ratio, fake
    printWriter.println(" 2");
    // membrane rois
    List<RefinementRoi> memRios = chomboSolverSpec.getMembraneRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Membrane + " " + memRios.size());
    for (RefinementRoi roi : memRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        // level tagsGrow ROIexpression
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    List<RefinementRoi> volRios = chomboSolverSpec.getVolumeRefinementRois();
    printWriter.println(FVInputFileKeyword.REFINEMENT_ROIS + " " + RoiType.Volume + " " + volRios.size());
    for (RefinementRoi roi : volRios) {
        if (roi.getRoiExpression() == null) {
            throw new SolverException("ROI expression cannot be null");
        }
        printWriter.println(roi.getLevel() + " " + roi.getRoiExpression().infix() + ";");
    }
    printWriter.println(FVInputFileKeyword.VIEW_LEVEL + " " + chomboSolverSpec.getViewLevel());
    printWriter.println(FVInputFileKeyword.CHOMBO_SPEC_END);
    printWriter.println();
}

Example 14 with Origin

use of org.vcell.util.Origin in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeDistanceMapFile.

private void writeDistanceMapFile(Vect3d deltaX, SubvolumeSignedDistanceMap distanceMap, File file) throws IOException {
    Origin origin = resampledGeometry.getOrigin();
    String[] headerNames = { "dimension", "Nx", "Ny", "Nz", "Dx", "Dy", "Dz", "firstX", "firstY", "firstZ" };
    double[] pointsX = distanceMap.getSamplesX();
    double[] pointsY = distanceMap.getSamplesY();
    double[] pointsZ = distanceMap.getSamplesZ();
    double[] distances = distanceMap.getSignedDistances();
    int Nx = pointsX.length;
    int Ny = pointsY.length;
    int Nz = pointsZ.length;
    double Dx = deltaX.getX();
    double Dy = deltaX.getY();
    double Dz = deltaX.getZ();
    assert Math.abs(origin.getX() - pointsX[0]) < 2e-8;
    assert Math.abs(origin.getY() - pointsY[0]) < 2e-8;
    // Chombo 2D uses middle slice, so skip assert in that case
    assert Math.abs(origin.getZ() - pointsZ[0]) < 2e-8 || bChomboSolver && resampledGeometry.getDimension() == 2 : "z not starting at origin";
    double[] header = new double[headerNames.length];
    int count = -1;
    header[++count] = resampledGeometry.getDimension();
    header[++count] = Nx;
    header[++count] = Ny;
    header[++count] = Nz;
    header[++count] = Dx;
    header[++count] = Dy;
    header[++count] = Dz;
    header[++count] = origin.getX();
    header[++count] = origin.getY();
    header[++count] = origin.getZ();
    DataOutputStream dos = null;
    try {
        dos = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(file)));
        for (double d : header) {
            dos.writeDouble(d);
        }
        for (double d : distances) {
            dos.writeDouble(d);
        }
    } finally {
        if (dos != null) {
            dos.close();
        }
    }
}

Example 15 with Origin

use of org.vcell.util.Origin in project vcell by virtualcell.

the class RunRefSimulationFastOp method runRefSimulation.

private RowColumnResultSet runRefSimulation(ROI cellROI, ROI[] imageDataROIs, UShortImage psf, FloatImage initRefConc, double experimentalRecoveryTime, LocalWorkspace localWorkspace, ClientTaskStatusSupport progressListener) throws Exception {
    User owner = LocalWorkspace.getDefaultOwner();
    KeyValue simKey = LocalWorkspace.createNewKeyValue();
    // 
    // save first image from normalized time series as the initial concentration field data
    // 
    ExternalDataInfo initialConcentrationExtData = createNewExternalDataInfo(localWorkspace, INITCONC_EXTDATA_NAME);
    Extent extent = initRefConc.getExtent();
    Origin origin = initRefConc.getOrigin();
    ISize isize = new ISize(initRefConc.getNumX(), initRefConc.getNumY(), initRefConc.getNumZ());
    saveExternalData(initRefConc, INITCONC_EXTDATA_VARNAME, initialConcentrationExtData.getExternalDataIdentifier(), localWorkspace);
    FieldFunctionArguments initConditionFFA = new FieldFunctionArguments(INITCONC_EXTDATA_NAME, INITCONC_EXTDATA_VARNAME, new Expression(0.0), VariableType.VOLUME);
    // 
    // save ROIs as a multivariate field data
    // 
    ExternalDataInfo roiExtData = createNewExternalDataInfo(localWorkspace, ROI_EXTDATA_NAME);
    saveROIsAsExternalData(imageDataROIs, localWorkspace, roiExtData.getExternalDataIdentifier());
    ArrayList<FieldFunctionArguments> roiFFAs = new ArrayList<FieldFunctionArguments>();
    for (ROI roi : imageDataROIs) {
        roiFFAs.add(new FieldFunctionArguments(ROI_EXTDATA_NAME, ROI_MASK_NAME_PREFIX + roi.getROIName(), new Expression(0.0), VariableType.VOLUME));
    }
    // 
    // save PSF as a field data
    // 
    ExternalDataInfo psfExtData = createNewExternalDataInfo(localWorkspace, PSF_EXTDATA_NAME);
    savePsfAsExternalData(psf, PSF_EXTDATA_VARNAME, psfExtData.getExternalDataIdentifier(), localWorkspace);
    FieldFunctionArguments psfFFA = new FieldFunctionArguments(PSF_EXTDATA_NAME, PSF_EXTDATA_VARNAME, new Expression(0.0), VariableType.VOLUME);
    TimeBounds timeBounds = getEstimatedRefTimeBound(experimentalRecoveryTime);
    double timeStepVal = REFERENCE_DIFF_DELTAT;
    Expression chirpedDiffusionRate = new Expression(REFERENCE_DIFF_RATE_COEFFICIENT + "*(t+" + REFERENCE_DIFF_DELTAT + ")");
    BioModel bioModel = createRefSimBioModel(simKey, owner, origin, extent, cellROI, timeStepVal, timeBounds, VAR_NAME, new Expression(initConditionFFA.infix()), psfFFA, chirpedDiffusionRate);
    if (progressListener != null) {
        progressListener.setMessage("Running Reference Simulation...");
    }
    // run simulation
    Simulation simulation = bioModel.getSimulation(0);
    ROIDataGenerator roiDataGenerator = getROIDataGenerator(localWorkspace, imageDataROIs);
    simulation.getMathDescription().getPostProcessingBlock().addDataGenerator(roiDataGenerator);
    runFVSolverStandalone(new File(localWorkspace.getDefaultSimDataDirectory()), simulation, initialConcentrationExtData.getExternalDataIdentifier(), roiExtData.getExternalDataIdentifier(), psfExtData.getExternalDataIdentifier(), progressListener, true);
    KeyValue referenceSimKeyValue = simulation.getVersion().getVersionKey();
    VCSimulationIdentifier vcSimID = new VCSimulationIdentifier(referenceSimKeyValue, LocalWorkspace.getDefaultOwner());
    VCSimulationDataIdentifier vcSimDataID = new VCSimulationDataIdentifier(vcSimID, 0);
    File hdf5File = new File(localWorkspace.getDefaultSimDataDirectory(), vcSimDataID.getID() + SimDataConstants.DATA_PROCESSING_OUTPUT_EXTENSION_HDF5);
    // get post processing info (time points, variable sizes)
    DataOperation.DataProcessingOutputInfoOP dataOperationInfo = new DataOperation.DataProcessingOutputInfoOP(null, /*no vcDataIdentifier OK*/
    false, null);
    DataOperationResults.DataProcessingOutputInfo dataProcessingOutputInfo = (DataOperationResults.DataProcessingOutputInfo) DataSetControllerImpl.getDataProcessingOutput(dataOperationInfo, hdf5File);
    // get post processing data
    DataOperation.DataProcessingOutputDataValuesOP dataOperationDataValues = new DataOperation.DataProcessingOutputDataValuesOP(null, /*no vcDataIdentifier OK*/
    ROI_EXTDATA_NAME, TimePointHelper.createAllTimeTimePointHelper(), DataIndexHelper.createSliceDataIndexHelper(0), null, null);
    DataOperationResults.DataProcessingOutputDataValues dataProcessingOutputDataValues = (DataOperationResults.DataProcessingOutputDataValues) DataSetControllerImpl.getDataProcessingOutput(dataOperationDataValues, hdf5File);
    // 
    // delete the simulation files
    // 
    // 
    // remove reference simulation files and field data files
    // 
    File userDir = new File(localWorkspace.getDefaultSimDataDirectory());
    File[] oldSimFilesToDelete = getSimulationFileNames(userDir, referenceSimKeyValue);
    for (int i = 0; oldSimFilesToDelete != null && i < oldSimFilesToDelete.length; i++) {
        oldSimFilesToDelete[i].delete();
    }
    deleteCanonicalExternalData(localWorkspace, initialConcentrationExtData.getExternalDataIdentifier());
    deleteCanonicalExternalData(localWorkspace, roiExtData.getExternalDataIdentifier());
    deleteCanonicalExternalData(localWorkspace, roiExtData.getExternalDataIdentifier());
    // get ref sim time points (distorted by time dilation acceleration)
    double[] rawRefDataTimePoints = dataProcessingOutputInfo.getVariableTimePoints();
    // get shifted time points
    double[] correctedRefDataTimePoints = shiftTimeForBaseDiffRate(rawRefDataTimePoints);
    double[][] refData = dataProcessingOutputDataValues.getDataValues();
    // 
    // for rowColumnResultSet with { "t", "roi1", .... , "roiN" } for reference data
    // 
    int numROIs = imageDataROIs.length;
    String[] columnNames = new String[numROIs + 1];
    columnNames[0] = "t";
    for (int i = 0; i < numROIs; i++) {
        columnNames[i + 1] = imageDataROIs[i].getROIName();
    }
    RowColumnResultSet reducedData = new RowColumnResultSet(columnNames);
    for (int i = 0; i < correctedRefDataTimePoints.length; i++) {
        double[] row = new double[numROIs + 1];
        row[0] = correctedRefDataTimePoints[i];
        double[] roiData = refData[i];
        for (int j = 0; j < numROIs; j++) {
            // roiData[0] is the average over the cell .. postbleach this shouldn't change for pure diffusion
            row[j + 1] = roiData[j + 1];
        }
        reducedData.addRow(row);
    }
    return reducedData;
}