Search in sources :

Example 51 with Origin

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

the class ITextWriter method writeGeom.

// Section used can be a chapter or a section of one, based on the document type.
protected void writeGeom(Section container, Geometry geom, GeometryContext geomCont) throws Exception {
    try {
        Section geomSection = container.addSection("Geometry: " + geom.getName(), container.numberDepth() + 1);
        if (geom.getDimension() == 0) {
            Paragraph p = new Paragraph(new Phrase("Non spatial geometry."));
            p.setAlignment(Paragraph.ALIGN_CENTER);
            geomSection.add(p);
            return;
        }
        ByteArrayOutputStream bos = generateGeometryImage(geom);
        com.lowagie.text.Image geomImage = com.lowagie.text.Image.getInstance(java.awt.Toolkit.getDefaultToolkit().createImage(bos.toByteArray()), null);
        geomImage.setAlignment(Table.ALIGN_LEFT);
        geomSection.add(geomImage);
        // addImage(geomSection, bos);
        Table geomTable = getTable(2, 50, 1, 2, 2);
        geomTable.setAlignment(Table.ALIGN_LEFT);
        Extent extent = geom.getExtent();
        String extentStr = "(" + extent.getX() + ", " + extent.getY() + ", " + extent.getZ() + ")";
        Origin origin = geom.getOrigin();
        String originStr = "(" + origin.getX() + ", " + origin.getY() + ", " + origin.getZ() + ")";
        geomTable.addCell(createCell("Size", getFont()));
        geomTable.addCell(createCell(extentStr, getFont()));
        geomTable.addCell(createCell("Origin", getFont()));
        geomTable.addCell(createCell(originStr, getFont()));
        geomSection.add(geomTable);
    } catch (Exception e) {
        System.err.println("Unable to add geometry image to report.");
        e.printStackTrace();
    }
}
Also used : Origin(org.vcell.util.Origin) Table(com.lowagie.text.Table) Extent(org.vcell.util.Extent) Phrase(com.lowagie.text.Phrase) ByteArrayOutputStream(java.io.ByteArrayOutputStream) SimulationContext(cbit.vcell.mapping.SimulationContext) GeometryContext(cbit.vcell.mapping.GeometryContext) ReactionContext(cbit.vcell.mapping.ReactionContext) Section(com.lowagie.text.Section) DocumentException(com.lowagie.text.DocumentException) ExpressionException(cbit.vcell.parser.ExpressionException) Paragraph(com.lowagie.text.Paragraph)

Example 52 with Origin

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

the class SBMLExporter method addGeometry.

private void addGeometry() throws SbmlException {
    SpatialModelPlugin mplugin = (SpatialModelPlugin) sbmlModel.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
    // Creates a geometry object via SpatialModelPlugin object.
    org.sbml.jsbml.ext.spatial.Geometry sbmlGeometry = mplugin.createGeometry();
    sbmlGeometry.setCoordinateSystem(GeometryKind.cartesian);
    sbmlGeometry.setSpatialId("vcell");
    Geometry vcGeometry = getSelectedSimContext().getGeometry();
    Model vcModel = getSelectedSimContext().getModel();
    // 
    // list of CoordinateComponents : 1 if geometry is 1-d, 2 if geometry is 2-d, 3 if geometry is 3-d
    // 
    int dimension = vcGeometry.getDimension();
    Extent vcExtent = vcGeometry.getExtent();
    Origin vcOrigin = vcGeometry.getOrigin();
    // add x coordinate component
    CoordinateComponent xComp = sbmlGeometry.createCoordinateComponent();
    xComp.setSpatialId(vcModel.getX().getName());
    xComp.setType(CoordinateKind.cartesianX);
    final UnitDefinition sbmlUnitDef_length = getOrCreateSBMLUnit(vcModel.getUnitSystem().getLengthUnit());
    xComp.setUnits(sbmlUnitDef_length);
    Boundary minX = new Boundary();
    xComp.setBoundaryMinimum(minX);
    minX.setSpatialId("Xmin");
    minX.setValue(vcOrigin.getX());
    Boundary maxX = new Boundary();
    xComp.setBoundaryMaximum(maxX);
    maxX.setSpatialId("Xmax");
    maxX.setValue(vcOrigin.getX() + (vcExtent.getX()));
    org.sbml.jsbml.Parameter pX = sbmlModel.createParameter();
    pX.setId(vcModel.getX().getName());
    pX.setValue(0.0);
    pX.setConstant(false);
    pX.setUnits(sbmlUnitDef_length);
    SpatialParameterPlugin spPluginPx = (SpatialParameterPlugin) pX.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
    SpatialSymbolReference spSymRefPx = new SpatialSymbolReference();
    spPluginPx.setParamType(spSymRefPx);
    spSymRefPx.setSpatialRef(xComp.getSpatialId());
    // add y coordinate component
    if (dimension == 2 || dimension == 3) {
        CoordinateComponent yComp = sbmlGeometry.createCoordinateComponent();
        yComp.setSpatialId(vcModel.getY().getName());
        yComp.setType(CoordinateKind.cartesianY);
        yComp.setUnits(sbmlUnitDef_length);
        Boundary minY = new Boundary();
        yComp.setBoundaryMinimum(minY);
        minY.setSpatialId("Ymin");
        minY.setValue(vcOrigin.getY());
        Boundary maxY = new Boundary();
        yComp.setBoundaryMaximum(maxY);
        maxY.setSpatialId("Ymax");
        maxY.setValue(vcOrigin.getY() + (vcExtent.getY()));
        org.sbml.jsbml.Parameter pY = sbmlModel.createParameter();
        pY.setId(vcModel.getY().getName());
        pY.setValue(0.0);
        pY.setConstant(false);
        pY.setUnits(sbmlUnitDef_length);
        SpatialParameterPlugin spPluginPy = (SpatialParameterPlugin) pY.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        SpatialSymbolReference spSymRefPy = new SpatialSymbolReference();
        spPluginPy.setParamType(spSymRefPy);
        spSymRefPy.setSpatialRef(yComp.getSpatialId());
    }
    // add z coordinate component
    if (dimension == 3) {
        CoordinateComponent zComp = sbmlGeometry.createCoordinateComponent();
        zComp.setSpatialId(vcModel.getZ().getName());
        zComp.setType(CoordinateKind.cartesianZ);
        zComp.setUnits(sbmlUnitDef_length);
        Boundary minZ = new Boundary();
        zComp.setBoundaryMinimum(minZ);
        minZ.setSpatialId("Zmin");
        minZ.setValue(vcOrigin.getZ());
        Boundary maxZ = new Boundary();
        zComp.setBoundaryMaximum(maxZ);
        maxZ.setSpatialId("Zmax");
        maxZ.setValue(vcOrigin.getZ() + (vcExtent.getZ()));
        org.sbml.jsbml.Parameter pZ = sbmlModel.createParameter();
        pZ.setId(vcModel.getZ().getName());
        pZ.setValue(0.0);
        pZ.setConstant(false);
        pZ.setUnits(sbmlUnitDef_length);
        SpatialParameterPlugin spPluginPz = (SpatialParameterPlugin) pZ.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        SpatialSymbolReference spSymRefPz = new SpatialSymbolReference();
        spPluginPz.setParamType(spSymRefPz);
        spSymRefPz.setSpatialRef(zComp.getSpatialId());
    }
    // 
    // list of compartmentMappings : VC structureMappings
    // 
    GeometryContext vcGeoContext = getSelectedSimContext().getGeometryContext();
    StructureMapping[] vcStrucMappings = vcGeoContext.getStructureMappings();
    for (int i = 0; i < vcStrucMappings.length; i++) {
        StructureMapping vcStructMapping = vcStrucMappings[i];
        String structName = vcStructMapping.getStructure().getName();
        Compartment comp = sbmlModel.getCompartment(TokenMangler.mangleToSName(structName));
        SpatialCompartmentPlugin cplugin = (SpatialCompartmentPlugin) comp.getPlugin(SBMLUtils.SBML_SPATIAL_NS_PREFIX);
        GeometryClass gc = vcStructMapping.getGeometryClass();
        if (!goodPointer(gc, GeometryClass.class, structName)) {
            continue;
        }
        CompartmentMapping compMapping = new CompartmentMapping();
        cplugin.setCompartmentMapping(compMapping);
        String geomClassName = gc.getName();
        String id = TokenMangler.mangleToSName(geomClassName + structName);
        compMapping.setSpatialId(id);
        compMapping.setDomainType(TokenMangler.mangleToSName(DOMAIN_TYPE_PREFIX + geomClassName));
        try {
            StructureMappingParameter usp = vcStructMapping.getUnitSizeParameter();
            Expression e = usp.getExpression();
            if (goodPointer(e, Expression.class, id)) {
                compMapping.setUnitSize(e.evaluateConstant());
            }
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to create compartment mapping for structureMapping '" + compMapping.getId() + "' : " + e.getMessage());
        }
    }
    // 
    // list of domain types : subvolumes and surface classes from VC
    // 
    boolean bAnyAnalyticSubvolumes = false;
    boolean bAnyImageSubvolumes = false;
    boolean bAnyCSGSubvolumes = false;
    GeometryClass[] vcGeomClasses = vcGeometry.getGeometryClasses();
    int numSubVols = 0;
    for (int i = 0; i < vcGeomClasses.length; i++) {
        DomainType domainType = sbmlGeometry.createDomainType();
        domainType.setSpatialId(DOMAIN_TYPE_PREFIX + vcGeomClasses[i].getName());
        if (vcGeomClasses[i] instanceof SubVolume) {
            if (((SubVolume) vcGeomClasses[i]) instanceof AnalyticSubVolume) {
                bAnyAnalyticSubvolumes = true;
            } else if (((SubVolume) vcGeomClasses[i]) instanceof ImageSubVolume) {
                bAnyImageSubvolumes = true;
            } else if (((SubVolume) vcGeomClasses[i]) instanceof CSGObject) {
                bAnyCSGSubvolumes = true;
            }
            domainType.setSpatialDimensions(3);
            numSubVols++;
        } else if (vcGeomClasses[i] instanceof SurfaceClass) {
            domainType.setSpatialDimensions(2);
        }
    }
    // 
    // list of domains, adjacent domains : from VC geometricRegions
    // 
    GeometrySurfaceDescription vcGSD = vcGeometry.getGeometrySurfaceDescription();
    if (vcGSD.getRegionImage() == null) {
        try {
            vcGSD.updateAll();
        } catch (Exception e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Unable to generate region images for geometry");
        }
    }
    GeometricRegion[] vcGeometricRegions = vcGSD.getGeometricRegions();
    ISize sampleSize = vcGSD.getVolumeSampleSize();
    int numX = sampleSize.getX();
    int numY = sampleSize.getY();
    int numZ = sampleSize.getZ();
    double ox = vcOrigin.getX();
    double oy = vcOrigin.getY();
    double oz = vcOrigin.getZ();
    RegionInfo[] regionInfos = vcGSD.getRegionImage().getRegionInfos();
    for (int i = 0; i < vcGeometricRegions.length; i++) {
        // domains
        Domain domain = sbmlGeometry.createDomain();
        domain.setSpatialId(vcGeometricRegions[i].getName());
        if (vcGeometricRegions[i] instanceof VolumeGeometricRegion) {
            domain.setDomainType(DOMAIN_TYPE_PREFIX + ((VolumeGeometricRegion) vcGeometricRegions[i]).getSubVolume().getName());
            // 
            // get a list of interior points ... should probably use the distance map to find a point
            // furthest inside (or several points associated with the morphological skeleton).
            // 
            InteriorPoint interiorPt = domain.createInteriorPoint();
            int regionID = ((VolumeGeometricRegion) vcGeometricRegions[i]).getRegionID();
            boolean bFound = false;
            int regInfoIndx = 0;
            for (int j = 0; j < regionInfos.length; j++) {
                regInfoIndx = j;
                if (regionInfos[j].getRegionIndex() == regionID) {
                    int volIndx = 0;
                    for (int z = 0; z < numZ && !bFound; z++) {
                        for (int y = 0; y < numY && !bFound; y++) {
                            for (int x = 0; x < numX && !bFound; x++) {
                                if (regionInfos[j].isIndexInRegion(volIndx)) {
                                    bFound = true;
                                    double unit_z = (numZ > 1) ? ((double) z) / (numZ - 1) : 0.5;
                                    double coordZ = oz + vcExtent.getZ() * unit_z;
                                    double unit_y = (numY > 1) ? ((double) y) / (numY - 1) : 0.5;
                                    double coordY = oy + vcExtent.getY() * unit_y;
                                    double unit_x = (numX > 1) ? ((double) x) / (numX - 1) : 0.5;
                                    double coordX = ox + vcExtent.getX() * unit_x;
                                    interiorPt.setCoord1(coordX);
                                    interiorPt.setCoord2(coordY);
                                    interiorPt.setCoord3(coordZ);
                                }
                                volIndx++;
                            }
                        // end - for x
                        }
                    // end - for y
                    }
                // end - for z
                }
            // end if
            }
            // end for regionInfos
            if (!bFound) {
                throw new RuntimeException("Unable to find interior point for region '" + regionInfos[regInfoIndx].toString());
            }
        } else if (vcGeometricRegions[i] instanceof SurfaceGeometricRegion) {
            SurfaceGeometricRegion vcSurfaceGeomReg = (SurfaceGeometricRegion) vcGeometricRegions[i];
            GeometricRegion geomRegion0 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[0];
            GeometricRegion geomRegion1 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[1];
            SurfaceClass surfaceClass = vcGSD.getSurfaceClass(((VolumeGeometricRegion) geomRegion0).getSubVolume(), ((VolumeGeometricRegion) geomRegion1).getSubVolume());
            domain.setDomainType(DOMAIN_TYPE_PREFIX + surfaceClass.getName());
            // adjacent domains : 2 adjacent domain objects for each surfaceClass in VC.
            // adjacent domain 1
            GeometricRegion adjGeomRegion0 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[0];
            GeometricRegion adjGeomRegion1 = vcSurfaceGeomReg.getAdjacentGeometricRegions()[1];
            AdjacentDomains adjDomain = new AdjacentDomains();
            adjDomain.setSpatialId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + adjGeomRegion0.getName()));
            adjDomain.setDomain1(vcSurfaceGeomReg.getName());
            adjDomain.setDomain2(adjGeomRegion0.getName());
            sbmlGeometry.addAdjacentDomain(adjDomain);
            // adj domain 2
            adjDomain = new AdjacentDomains();
            adjDomain.setSpatialId(TokenMangler.mangleToSName(vcSurfaceGeomReg.getName() + "_" + adjGeomRegion1.getName()));
            adjDomain.setDomain1(vcSurfaceGeomReg.getName());
            adjDomain.setDomain2(adjGeomRegion1.getName());
            sbmlGeometry.addAdjacentDomain(adjDomain);
        }
    }
    // 
    if (bAnyAnalyticSubvolumes && !bAnyImageSubvolumes && !bAnyCSGSubvolumes) {
        AnalyticGeometry sbmlAnalyticGeomDefinition = sbmlGeometry.createAnalyticGeometry();
        sbmlAnalyticGeomDefinition.setSpatialId(TokenMangler.mangleToSName("Analytic_" + vcGeometry.getName()));
        sbmlAnalyticGeomDefinition.setIsActive(true);
        for (int i = 0; i < vcGeomClasses.length; i++) {
            if (vcGeomClasses[i] instanceof AnalyticSubVolume) {
                AnalyticVolume analyticVol = sbmlAnalyticGeomDefinition.createAnalyticVolume();
                analyticVol.setSpatialId(vcGeomClasses[i].getName());
                analyticVol.setDomainType(DOMAIN_TYPE_PREFIX + vcGeomClasses[i].getName());
                analyticVol.setFunctionType(FunctionKind.layered);
                analyticVol.setOrdinal(numSubVols - (i + 1));
                Expression expr = ((AnalyticSubVolume) vcGeomClasses[i]).getExpression();
                try {
                    String mathMLStr = ExpressionMathMLPrinter.getMathML(expr, true, MathType.BOOLEAN);
                    ASTNode mathMLNode = ASTNode.readMathMLFromString(mathMLStr);
                    analyticVol.setMath(mathMLNode);
                } catch (Exception e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Error converting VC subvolume expression to mathML" + e.getMessage());
                }
            }
        }
    }
    // 
    if (!bAnyAnalyticSubvolumes && !bAnyImageSubvolumes && bAnyCSGSubvolumes) {
        CSGeometry sbmlCSGeomDefinition = new CSGeometry();
        sbmlGeometry.addGeometryDefinition(sbmlCSGeomDefinition);
        sbmlCSGeomDefinition.setSpatialId(TokenMangler.mangleToSName("CSG_" + vcGeometry.getName()));
        for (int i = 0; i < vcGeomClasses.length; i++) {
            if (vcGeomClasses[i] instanceof CSGObject) {
                CSGObject vcellCSGObject = (CSGObject) vcGeomClasses[i];
                org.sbml.jsbml.ext.spatial.CSGObject sbmlCSGObject = new org.sbml.jsbml.ext.spatial.CSGObject();
                sbmlCSGeomDefinition.addCSGObject(sbmlCSGObject);
                sbmlCSGObject.setSpatialId(vcellCSGObject.getName());
                sbmlCSGObject.setDomainType(DOMAIN_TYPE_PREFIX + vcellCSGObject.getName());
                // the ordinal should the the least for the default/background subVolume
                sbmlCSGObject.setOrdinal(numSubVols - (i + 1));
                org.sbml.jsbml.ext.spatial.CSGNode sbmlcsgNode = getSBMLCSGNode(vcellCSGObject.getRoot());
                sbmlCSGObject.setCSGNode(sbmlcsgNode);
            }
        }
    }
    // 
    // add "Segmented" and "DistanceMap" SampledField Geometries
    // 
    final boolean bVCGeometryIsImage = bAnyImageSubvolumes && !bAnyAnalyticSubvolumes && !bAnyCSGSubvolumes;
    // 55if (bAnyAnalyticSubvolumes || bAnyImageSubvolumes || bAnyCSGSubvolumes){
    if (bVCGeometryIsImage) {
        // 
        // add "Segmented" SampledFieldGeometry
        // 
        SampledFieldGeometry segmentedImageSampledFieldGeometry = sbmlGeometry.createSampledFieldGeometry();
        segmentedImageSampledFieldGeometry.setSpatialId(TokenMangler.mangleToSName("SegmentedImage_" + vcGeometry.getName()));
        segmentedImageSampledFieldGeometry.setIsActive(true);
        // 55boolean bVCGeometryIsImage = bAnyImageSubvolumes && !bAnyAnalyticSubvolumes && !bAnyCSGSubvolumes;
        Geometry vcImageGeometry = null;
        {
            if (bVCGeometryIsImage) {
                // make a resampled image;
                if (dimension == 3) {
                    try {
                        ISize imageSize = vcGeometry.getGeometrySpec().getDefaultSampledImageSize();
                        vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT());
                        vcImageGeometry = RayCaster.resampleGeometry(new GeometryThumbnailImageFactoryAWT(), vcGeometry, imageSize);
                    } catch (Throwable e) {
                        e.printStackTrace(System.out);
                        throw new RuntimeException("Unable to convert the original analytic or constructed solid geometry to image-based geometry : " + e.getMessage());
                    }
                } else {
                    try {
                        vcGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, false);
                        GeometrySpec origGeometrySpec = vcGeometry.getGeometrySpec();
                        VCImage newVCImage = origGeometrySpec.getSampledImage().getCurrentValue();
                        // 
                        // construct the new geometry with the sampled VCImage.
                        // 
                        vcImageGeometry = new Geometry(vcGeometry.getName() + "_asImage", newVCImage);
                        vcImageGeometry.getGeometrySpec().setExtent(vcGeometry.getExtent());
                        vcImageGeometry.getGeometrySpec().setOrigin(vcGeometry.getOrigin());
                        vcImageGeometry.setDescription(vcGeometry.getDescription());
                        vcImageGeometry.getGeometrySurfaceDescription().setFilterCutoffFrequency(vcGeometry.getGeometrySurfaceDescription().getFilterCutoffFrequency());
                        vcImageGeometry.precomputeAll(new GeometryThumbnailImageFactoryAWT(), true, true);
                    } catch (Exception e) {
                        e.printStackTrace(System.out);
                        throw new RuntimeException("Unable to convert the original analytic or constructed solid geometry to image-based geometry : " + e.getMessage());
                    }
                }
                GeometryClass[] vcImageGeomClasses = vcImageGeometry.getGeometryClasses();
                for (int j = 0; j < vcImageGeomClasses.length; j++) {
                    if (vcImageGeomClasses[j] instanceof ImageSubVolume) {
                        SampledVolume sampledVol = segmentedImageSampledFieldGeometry.createSampledVolume();
                        sampledVol.setSpatialId(vcGeomClasses[j].getName());
                        sampledVol.setDomainType(DOMAIN_TYPE_PREFIX + vcGeomClasses[j].getName());
                        sampledVol.setSampledValue(((ImageSubVolume) vcImageGeomClasses[j]).getPixelValue());
                    }
                }
                // add sampledField to sampledFieldGeometry
                SampledField segmentedImageSampledField = sbmlGeometry.createSampledField();
                VCImage vcImage = vcImageGeometry.getGeometrySpec().getImage();
                segmentedImageSampledField.setSpatialId("SegmentedImageSampledField");
                segmentedImageSampledField.setNumSamples1(vcImage.getNumX());
                segmentedImageSampledField.setNumSamples2(vcImage.getNumY());
                segmentedImageSampledField.setNumSamples3(vcImage.getNumZ());
                segmentedImageSampledField.setInterpolationType(InterpolationKind.nearestneighbor);
                segmentedImageSampledField.setCompression(CompressionKind.uncompressed);
                segmentedImageSampledField.setDataType(DataKind.UINT8);
                segmentedImageSampledFieldGeometry.setSampledField(segmentedImageSampledField.getId());
                try {
                    byte[] vcImagePixelsBytes = vcImage.getPixels();
                    // imageData.setCompression("");
                    StringBuffer sb = new StringBuffer();
                    for (int i = 0; i < vcImagePixelsBytes.length; i++) {
                        int uint8_sample = ((int) vcImagePixelsBytes[i]) & 0xff;
                        sb.append(uint8_sample + " ");
                    }
                    segmentedImageSampledField.setSamplesLength(vcImage.getNumXYZ());
                    segmentedImageSampledField.setSamples(sb.toString().trim());
                } catch (ImageException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Unable to export image from VCell to SBML : " + e.getMessage());
                }
            }
        }
    /*		
		//
		// add "DistanceMap" SampledFieldGeometry if there are exactly two subvolumes (else need more fields) and geometry is 3d.
		//
		if (numSubVols==2 && dimension == 3){
			SignedDistanceMap[] distanceMaps = null;
			try {
				distanceMaps = DistanceMapGenerator.computeDistanceMaps(vcImageGeometry, vcImageGeometry.getGeometrySpec().getImage(), false, false);
			} catch (ImageException e) {
				e.printStackTrace(System.out);
				System.err.println("Unable to export distance map sampled field from VCell to SBML : " + e.getMessage());
				// throw new RuntimeException("Unable to export distance map sampled field from VCell to SBML : " + e.getMessage());
				
				// don't want to throw an exception and stop export because distance map geometry couldn't be exported. 
				// just 'return' from method (since this is the last thing that is being done in this method).
				return;
			}
			//
			// the two distanceMaps should be redundant (one is negation of the other) ... so choose first one for field.
			//
			double[] signedDistances = distanceMaps[0].getSignedDistances();
			SampledFieldGeometry distanceMapSampledFieldGeometry = sbmlGeometry.createSampledFieldGeometry();
			distanceMapSampledFieldGeometry.setSpatialId(TokenMangler.mangleToSName("DistanceMap_"+vcGeometry.getName()));
			SampledField distanceMapSampledField = distanceMapSampledFieldGeometry.createSampledField();
			distanceMapSampledField.setSpatialId("DistanceMapSampledField");
			distanceMapSampledField.setNumSamples1(distanceMaps[0].getSamplesX().length);
			distanceMapSampledField.setNumSamples2(distanceMaps[0].getSamplesY().length);
			distanceMapSampledField.setNumSamples3(distanceMaps[0].getSamplesZ().length);
			distanceMapSampledField.setDataType("real");
System.err.println("do we need distanceMapSampleField.setDataType()?");
			distanceMapSampledField.setInterpolationType("linear");
			ImageData distanceMapImageData = distanceMapSampledField.createImageData();
			distanceMapImageData.setDataType("int16");
System.err.println("should be:\n  distanceMapImageData.setDataType(\"float32\")");
//					distanceMapImageData.setCompression("");

			double maxAbsValue = 0;
			for (int i = 0; i < signedDistances.length; i++) {
				maxAbsValue = Math.max(maxAbsValue,Math.abs(signedDistances[i]));
			}
			if (maxAbsValue==0.0){
				throw new RuntimeException("computed distance map all zeros");
			}
			double scale = (Short.MAX_VALUE-1)/maxAbsValue;
			int[] scaledIntegerDistanceMap = new int[signedDistances.length];
			for (int i = 0; i < signedDistances.length; i++) {
				scaledIntegerDistanceMap[i] = (int)(scale * signedDistances[i]);
			}
			distanceMapImageData.setSamples(scaledIntegerDistanceMap, signedDistances.length);
System.err.println("should be:\n  distanceMapImageData.setSamples((float[])signedDistances,signedDistances.length)");
			SampledVolume sampledVol = distanceMapSampledFieldGeometry.createSampledVolume();
			sampledVol.setSpatialId(distanceMaps[0].getInsideSubvolumeName());
			sampledVol.setDomainType(DOMAIN_TYPE_PREFIX+distanceMaps[0].getInsideSubvolumeName());
			sampledVol.setSampledValue(255);
			sampledVol = distanceMapSampledFieldGeometry.createSampledVolume();
			sampledVol.setSpatialId(distanceMaps[1].getInsideSubvolumeName());
			sampledVol.setDomainType(DOMAIN_TYPE_PREFIX+distanceMaps[1].getInsideSubvolumeName());
			sampledVol.setSampledValue(1);
		}
*/
    }
// 
// add "SurfaceMesh" ParametricGeometry
// 
// if (bAnyAnalyticSubvolumes || bAnyImageSubvolumes || bAnyCSGSubvolumes){
// ParametricGeometry sbmlParametricGeomDefinition = sbmlGeometry.createParametricGeometry();
// sbmlParametricGeomDefinition.setSpatialId(TokenMangler.mangleToSName("SurfaceMesh_"+vcGeometry.getName()));
// xxxx
// }
}
Also used : Origin(org.vcell.util.Origin) CompartmentMapping(org.sbml.jsbml.ext.spatial.CompartmentMapping) Compartment(org.sbml.jsbml.Compartment) SpatialParameterPlugin(org.sbml.jsbml.ext.spatial.SpatialParameterPlugin) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) ExpressionException(cbit.vcell.parser.ExpressionException) Boundary(org.sbml.jsbml.ext.spatial.Boundary) GeometrySpec(cbit.vcell.geometry.GeometrySpec) DomainType(org.sbml.jsbml.ext.spatial.DomainType) SampledVolume(org.sbml.jsbml.ext.spatial.SampledVolume) SubVolume(cbit.vcell.geometry.SubVolume) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) GeometryContext(cbit.vcell.mapping.GeometryContext) SpatialCompartmentPlugin(org.sbml.jsbml.ext.spatial.SpatialCompartmentPlugin) CoordinateComponent(org.sbml.jsbml.ext.spatial.CoordinateComponent) SimulationContext(cbit.vcell.mapping.SimulationContext) SpeciesContext(cbit.vcell.model.SpeciesContext) GeometryContext(cbit.vcell.mapping.GeometryContext) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion) SpatialSymbolReference(org.sbml.jsbml.ext.spatial.SpatialSymbolReference) AnalyticVolume(org.sbml.jsbml.ext.spatial.AnalyticVolume) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) SampledField(org.sbml.jsbml.ext.spatial.SampledField) Domain(org.sbml.jsbml.ext.spatial.Domain) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) GeometryClass(cbit.vcell.geometry.GeometryClass) ImageException(cbit.image.ImageException) GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription) Extent(org.vcell.util.Extent) SurfaceClass(cbit.vcell.geometry.SurfaceClass) ISize(org.vcell.util.ISize) CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry) RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo) VCImage(cbit.image.VCImage) StructureMapping(cbit.vcell.mapping.StructureMapping) GeometryThumbnailImageFactoryAWT(cbit.vcell.geometry.GeometryThumbnailImageFactoryAWT) ASTNode(org.sbml.jsbml.ASTNode) CSGObject(cbit.vcell.geometry.CSGObject) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) UnitDefinition(org.sbml.jsbml.UnitDefinition) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) SpatialModelPlugin(org.sbml.jsbml.ext.spatial.SpatialModelPlugin) InteriorPoint(org.sbml.jsbml.ext.spatial.InteriorPoint) XMLStreamException(javax.xml.stream.XMLStreamException) SbmlException(org.vcell.sbml.SbmlException) ImageException(cbit.image.ImageException) SBMLException(org.sbml.jsbml.SBMLException) ExpressionException(cbit.vcell.parser.ExpressionException) AdjacentDomains(org.sbml.jsbml.ext.spatial.AdjacentDomains) Geometry(cbit.vcell.geometry.Geometry) SampledFieldGeometry(org.sbml.jsbml.ext.spatial.SampledFieldGeometry) AnalyticGeometry(org.sbml.jsbml.ext.spatial.AnalyticGeometry) CSGeometry(org.sbml.jsbml.ext.spatial.CSGeometry) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) BioModel(cbit.vcell.biomodel.BioModel)

Example 53 with Origin

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

the class SmoldynFileWriter method writeInitialConcentration.

private int writeInitialConcentration(ParticleInitialConditionConcentration initialConcentration, SubDomain subDomain, Variable variable, String variableName, StringBuilder sb) throws ExpressionException, MathException {
    SimpleSymbolTable simpleSymbolTable = new SimpleSymbolTable(new String[] { ReservedVariable.X.getName(), ReservedVariable.Y.getName(), ReservedVariable.Z.getName() });
    Expression disExpression = new Expression(initialConcentration.getDistribution());
    disExpression.bindExpression(simulationSymbolTable);
    disExpression = simulationSymbolTable.substituteFunctions(disExpression).flatten();
    disExpression.bindExpression(simpleSymbolTable);
    double[] values = new double[3];
    if (dimension == 1) {
        if (disExpression.getSymbolBinding(ReservedVariable.Y.getName()) != null || disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'y' or 'z'", dimension, variable, disExpression));
        }
    } else if (dimension == 2) {
        if (disExpression.getSymbolBinding(ReservedVariable.Z.getName()) != null) {
            throw new MathException(VCellErrorMessages.getSmoldynWrongCoordinates("'z'", dimension, variable, disExpression));
        }
    }
    int totalCount = 0;
    StringBuilder localsb = new StringBuilder();
    if (subDomain instanceof CompartmentSubDomain) {
        MeshSpecification meshSpecification = simulation.getMeshSpecification();
        ISize sampleSize = meshSpecification.getSamplingSize();
        int numX = sampleSize.getX();
        int numY = dimension < 2 ? 1 : sampleSize.getY();
        int numZ = dimension < 3 ? 1 : sampleSize.getZ();
        boolean bCellCentered = simulation.hasCellCenteredMesh();
        double dx = meshSpecification.getDx(bCellCentered);
        double dy = meshSpecification.getDy(bCellCentered);
        double dz = meshSpecification.getDz(bCellCentered);
        Origin origin = resampledGeometry.getGeometrySpec().getOrigin();
        double ox = origin.getX();
        double oy = origin.getY();
        double oz = origin.getZ();
        Extent extent = resampledGeometry.getExtent();
        double ex = extent.getX();
        double ey = extent.getY();
        double ez = extent.getZ();
        int offset = 0;
        for (int k = 0; k < numZ; k++) {
            double centerz = oz + k * dz;
            double loz = Math.max(oz, centerz - dz / 2);
            double hiz = Math.min(oz + ez, centerz + dz / 2);
            double lz = hiz - loz;
            values[2] = centerz;
            for (int j = 0; j < numY; j++) {
                double centery = oy + j * dy;
                double loy = Math.max(oy, centery - dy / 2);
                double hiy = Math.min(oy + ey, centery + dy / 2);
                values[1] = centery;
                double ly = hiy - loy;
                for (int i = 0; i < numX; i++) {
                    int regionIndex = resampledGeometry.getGeometrySurfaceDescription().getRegionImage().getRegionInfoFromOffset(offset).getRegionIndex();
                    offset++;
                    GeometricRegion region = resampledGeometry.getGeometrySurfaceDescription().getGeometricRegions(regionIndex);
                    if (region instanceof VolumeGeometricRegion) {
                        if (!((VolumeGeometricRegion) region).getSubVolume().getName().equals(subDomain.getName())) {
                            continue;
                        }
                    }
                    double centerx = ox + i * dx;
                    double lox = Math.max(ox, centerx - dx / 2);
                    double hix = Math.min(ox + ex, centerx + dx / 2);
                    double lx = hix - lox;
                    values[0] = centerx;
                    double volume = lx;
                    if (dimension > 1) {
                        volume *= ly;
                        if (dimension > 2) {
                            volume *= lz;
                        }
                    }
                    double expectedCount = disExpression.evaluateVector(values) * volume;
                    if (expectedCount <= 0) {
                        continue;
                    }
                    long count = dist.nextPoisson(expectedCount);
                    if (count <= 0) {
                        continue;
                    }
                    totalCount += count;
                    localsb.append(SmoldynVCellMapper.SmoldynKeyword.mol + " " + count + " " + variableName + " " + (float) lox + "-" + (float) hix);
                    if (lg.isDebugEnabled()) {
                        lg.debug("Component subdomain " + variableName + " count " + count);
                    }
                    if (dimension > 1) {
                        localsb.append(" " + loy + "-" + hiy);
                        if (dimension > 2) {
                            localsb.append(" " + loz + "-" + hiz);
                        }
                    }
                    localsb.append("\n");
                }
            }
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.compartment_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    } else if (subDomain instanceof MembraneSubDomain) {
        ArrayList<TrianglePanel> trianglePanelList = membraneSubdomainTriangleMap.get(subDomain);
        for (TrianglePanel trianglePanel : trianglePanelList) {
            Triangle triangle = trianglePanel.triangle;
            switch(dimension) {
                case 1:
                    values[0] = triangle.getNodes(0).getX();
                    break;
                case 2:
                    {
                        double centroidX = triangle.getNodes(0).getX();
                        double centroidY = triangle.getNodes(0).getY();
                        if (triangle.getNodes(0).getX() == triangle.getNodes(1).getX() && triangle.getNodes(0).getY() == triangle.getNodes(1).getY()) {
                            centroidX += triangle.getNodes(2).getX();
                            centroidY += triangle.getNodes(2).getY();
                        } else {
                            centroidX += triangle.getNodes(1).getX();
                            centroidY += triangle.getNodes(1).getY();
                        }
                        values[0] = centroidX / 2;
                        values[1] = centroidY / 2;
                        break;
                    }
                case 3:
                    {
                        double centroidX = triangle.getNodes(0).getX() + triangle.getNodes(1).getX() + triangle.getNodes(2).getX();
                        double centroidY = triangle.getNodes(0).getY() + triangle.getNodes(1).getY() + triangle.getNodes(2).getY();
                        double centroidZ = triangle.getNodes(0).getZ() + triangle.getNodes(1).getZ() + triangle.getNodes(2).getZ();
                        values[0] = centroidX / 3;
                        values[1] = centroidY / 3;
                        values[2] = centroidZ / 3;
                        break;
                    }
            }
            double expectedCount = disExpression.evaluateVector(values) * triangle.getArea();
            if (expectedCount <= 0) {
                continue;
            }
            long count = dist.nextPoisson(expectedCount);
            if (count <= 0) {
                continue;
            }
            totalCount += count;
            if (lg.isDebugEnabled()) {
                lg.debug("Membrane subdomain " + subDomain.getName() + ' ' + variableName + " count " + count);
            }
            localsb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol + " " + count + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.tri + " " + trianglePanel.name + "\n");
        }
        // otherwise we append the distributed molecules in different small boxes
        try {
            subsituteFlattenToConstant(disExpression);
            sb.append(SmoldynVCellMapper.SmoldynKeyword.surface_mol);
            sb.append(" " + totalCount + " " + variableName + " " + subDomain.getName() + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.all + "\n");
        } catch (// can not be evaluated to a constant
        Exception e) {
            sb.append(localsb);
        }
    }
    if (lg.isDebugEnabled()) {
        lg.debug("Subdomain " + subDomain.getName() + ' ' + variableName + " total count " + totalCount);
    }
    return totalCount;
}
Also used : Origin(org.vcell.util.Origin) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Extent(org.vcell.util.Extent) ISize(org.vcell.util.ISize) ArrayList(java.util.ArrayList) Triangle(cbit.vcell.geometry.surface.Triangle) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) MeshSpecification(cbit.vcell.solver.MeshSpecification) ProgrammingException(org.vcell.util.ProgrammingException) GeometryException(cbit.vcell.geometry.GeometryException) IOException(java.io.IOException) DataAccessException(org.vcell.util.DataAccessException) PropertyVetoException(java.beans.PropertyVetoException) DivideByZeroException(cbit.vcell.parser.DivideByZeroException) ImageException(cbit.image.ImageException) ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException) SolverException(cbit.vcell.solver.SolverException) ExpressionException(cbit.vcell.parser.ExpressionException) MathException(cbit.vcell.math.MathException) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion) SimpleSymbolTable(cbit.vcell.parser.SimpleSymbolTable) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)

Example 54 with Origin

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

the class SmoldynFileWriter method writeWallSurfaces.

private void writeWallSurfaces() throws SolverException {
    GeometrySurfaceDescription geometrySurfaceDescription = resampledGeometry.getGeometrySurfaceDescription();
    GeometrySpec geometrySpec = resampledGeometry.getGeometrySpec();
    SubVolume[] subVolumes = geometrySpec.getSubVolumes();
    printWriter.println("# boundaries");
    Origin origin = geometrySpec.getOrigin();
    Extent extent = geometrySpec.getExtent();
    Coordinate lowWall = new Coordinate(origin.getX(), origin.getY(), origin.getZ());
    Coordinate highWall = new Coordinate(origin.getX() + extent.getX(), origin.getY() + extent.getY(), origin.getZ() + extent.getZ());
    // potential artifact.
    if (bHasNoSurface) {
        SubDomain subDomain0 = mathDesc.getSubDomains().nextElement();
        CompartmentSubDomain compartSubDomain0 = null;
        compartSubDomain0 = (CompartmentSubDomain) subDomain0;
        // x
        if (compartSubDomain0.getBoundaryConditionXm().isPERIODIC()) {
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
        } else {
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 0 " + lowWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 0 " + highWall.getX() + " " + (compartSubDomain0.getBoundaryConditionXp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
        }
        if (dimension > 1) {
            // y
            if (compartSubDomain0.getBoundaryConditionYm().isPERIODIC()) {
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
            } else {
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 1 " + lowWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 1 " + highWall.getY() + " " + (compartSubDomain0.getBoundaryConditionYp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
            }
            if (dimension > 2) {
                // z
                if (compartSubDomain0.getBoundaryConditionZm().isPERIODIC()) {
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ() + " " + SmoldynVCellMapper.SmoldynKeyword.p);
                } else {
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.low_wall + " 2 " + lowWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZm().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.high_wall + " 2 " + highWall.getZ() + " " + (compartSubDomain0.getBoundaryConditionZp().isNEUMANN() ? SmoldynVCellMapper.SmoldynKeyword.r : SmoldynVCellMapper.SmoldynKeyword.a));
                }
            }
        }
        printWriter.println();
    } else {
        // x
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 0 " + lowWall.getX() + " " + highWall.getX());
        if (dimension > 1) {
            // y
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 1 " + lowWall.getY() + " " + highWall.getY());
            if (dimension > 2) {
                // z
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.boundaries + " 2 " + lowWall.getZ() + " " + highWall.getZ());
            }
        }
        printWriter.println();
        // bounding walls as surfaces
        // have to find boundary condition type
        ISize sampleSize = simulation.getMeshSpecification().getSamplingSize();
        int numX = sampleSize.getX();
        int numY = dimension < 2 ? 1 : sampleSize.getY();
        int numZ = dimension < 3 ? 1 : sampleSize.getZ();
        if (dimension > 2) {
            int[] k_wall = new int[] { 0, numZ - 1 };
            for (int k = 0; k < k_wall.length; k++) {
                for (int j = 0; j < numY; j++) {
                    for (int i = 0; i < numX; i++) {
                        int volIndex = k_wall[k] * numX * numY + j * numX + i;
                        for (SubVolume sv : subVolumes) {
                            // gather all the points in all the regions
                            GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                            RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                            for (GeometricRegion gr : geometricRegions) {
                                VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                                for (RegionInfo ri : regionInfos) {
                                    if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                        boundaryZSubVolumes.add(sv);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        if (dimension > 1) {
            int[] j_wall = new int[] { 0, numY - 1 };
            for (int k = 0; k < numZ; k++) {
                for (int j = 0; j < j_wall.length; j++) {
                    for (int i = 0; i < numX; i++) {
                        int volIndex = k * numX * numY + j_wall[j] * numX + i;
                        for (SubVolume sv : subVolumes) {
                            // gather all the points in all the regions
                            GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                            RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                            for (GeometricRegion gr : geometricRegions) {
                                VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                                for (RegionInfo ri : regionInfos) {
                                    if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                        boundaryYSubVolumes.add(sv);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        int[] i_wall = new int[] { 0, numX - 1 };
        for (int k = 0; k < numZ; k++) {
            for (int j = 0; j < numY; j++) {
                for (int i = 0; i < i_wall.length; i++) {
                    int volIndex = k * numX * numY + j * numX + i_wall[i];
                    for (SubVolume sv : subVolumes) {
                        // gather all the points in all the regions
                        GeometricRegion[] geometricRegions = geometrySurfaceDescription.getGeometricRegions(sv);
                        RegionInfo[] regionInfos = geometrySurfaceDescription.getRegionImage().getRegionInfos();
                        for (GeometricRegion gr : geometricRegions) {
                            VolumeGeometricRegion vgr = (VolumeGeometricRegion) gr;
                            for (RegionInfo ri : regionInfos) {
                                if (ri.getRegionIndex() == vgr.getRegionID() && ri.isIndexInRegion(volIndex)) {
                                    boundaryXSubVolumes.add(sv);
                                }
                            }
                        }
                    }
                }
            }
        }
        Set<SubVolume> boundarySubVolumes = new HashSet<SubVolume>();
        boundarySubVolumes.addAll(boundaryXSubVolumes);
        boundarySubVolumes.addAll(boundaryYSubVolumes);
        boundarySubVolumes.addAll(boundaryZSubVolumes);
        BoundaryConditionType[] computedBct = new BoundaryConditionType[dimension * 2];
        String[] smoldynBct = new String[dimension * 2];
        String[] wallNames = new String[] { "Xm", "Xp", "Ym", "Yp", "Zm", "Zp" };
        if (boundarySubVolumes.size() >= 1) {
            for (SubVolume sv : boundarySubVolumes) {
                CompartmentSubDomain csd = (CompartmentSubDomain) mathDesc.getSubDomain(sv.getName());
                BoundaryConditionType[] bct = new BoundaryConditionType[] { csd.getBoundaryConditionXm(), csd.getBoundaryConditionXp(), csd.getBoundaryConditionYm(), csd.getBoundaryConditionYp(), csd.getBoundaryConditionZm(), csd.getBoundaryConditionZp() };
                if (computedBct[0] == null) {
                    System.arraycopy(bct, 0, computedBct, 0, dimension * 2);
                    for (int i = 0; i < dimension * 2; i++) {
                        if (computedBct[i].isPERIODIC()) {
                            throw new SolverException("Models with both surfaces and periodic boundary conditions are not supported yet.");
                        }
                        smoldynBct[i] = computedBct[i].isDIRICHLET() ? SmoldynVCellMapper.SmoldynKeyword.absorb.name() : SmoldynVCellMapper.SmoldynKeyword.reflect.name();
                    }
                } else {
                    for (int i = 0; i < dimension * 2; i++) {
                        if (!computedBct[i].compareEqual(bct[i])) {
                            throw new SolverException(wallNames[i] + " wall has different boundary conditions");
                        }
                    }
                }
            }
        }
        printWriter.println("# bounding wall surface");
        // X walls
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_X);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[0]);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[1]);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
        // yz walls
        switch(dimension) {
            case 1:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX());
                break;
            case 2:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + extent.getY());
                break;
            case 3:
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +0 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -0 " + highWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getY() + " " + extent.getZ());
                break;
        }
        printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
        printWriter.println();
        if (dimension > 1) {
            // Y walls
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Y);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[2]);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[3]);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
            // xz walls
            switch(dimension) {
                case 2:
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + extent.getX());
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + extent.getX());
                    break;
                case 3:
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +1 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -1 " + lowWall.getX() + " " + highWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getZ());
                    break;
            }
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
            printWriter.println();
            if (dimension > 2) {
                // Z walls
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.start_surface + " " + VCellSmoldynKeyword.bounding_wall_surface_Z);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + "(" + SmoldynVCellMapper.SmoldynKeyword.up + ") " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.reflect);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.front + " " + smoldynBct[4]);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.action + " " + SmoldynVCellMapper.SmoldynKeyword.all + " " + SmoldynVCellMapper.SmoldynKeyword.back + " " + smoldynBct[5]);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.color + " " + SmoldynVCellMapper.SmoldynKeyword.both + " 1 1 1");
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.polygon + " " + SmoldynVCellMapper.SmoldynKeyword.both + " " + SmoldynVCellMapper.SmoldynKeyword.edge);
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_panels + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " 2");
                // xy walls
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " +2 " + lowWall.getX() + " " + lowWall.getY() + " " + lowWall.getZ() + " " + extent.getX() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.panel + " " + SmoldynVCellMapper.SmoldynKeyword.rect + " -2 " + lowWall.getX() + " " + lowWall.getY() + " " + highWall.getZ() + " " + extent.getX() + " " + extent.getY());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.end_surface);
                printWriter.println();
            }
        }
    }
}
Also used : Origin(org.vcell.util.Origin) GeometrySurfaceDescription(cbit.vcell.geometry.surface.GeometrySurfaceDescription) Extent(org.vcell.util.Extent) ISize(org.vcell.util.ISize) RegionInfo(cbit.vcell.geometry.RegionImage.RegionInfo) BoundaryConditionType(cbit.vcell.math.BoundaryConditionType) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) SurfaceGeometricRegion(cbit.vcell.geometry.surface.SurfaceGeometricRegion) VolumeGeometricRegion(cbit.vcell.geometry.surface.VolumeGeometricRegion) GeometricRegion(cbit.vcell.geometry.surface.GeometricRegion) GeometrySpec(cbit.vcell.geometry.GeometrySpec) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Coordinate(org.vcell.util.Coordinate) SubVolume(cbit.vcell.geometry.SubVolume) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SolverException(cbit.vcell.solver.SolverException) HashSet(java.util.HashSet)

Example 55 with Origin

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

the class RunRefSimulationFastOp method getROIDataGenerator.

private ROIDataGenerator getROIDataGenerator(LocalContext localWorkspace, ROI[] rois) throws ImageException, IOException {
    // create ROI image
    short[] roiFieldData = null;
    if (rois.length > 0) {
        Origin origin = new Origin(0, 0, 0);
        Extent extent = rois[0].getRoiImages()[0].getExtent();
        ISize isize = rois[0].getISize();
        int numROIX = rois[0].getISize().getX();
        int numROIY = rois[0].getISize().getY();
        roiFieldData = new short[numROIX * numROIY];
        short regionCounter = 1;
        for (int roiIdx = 0; roiIdx < rois.length; roiIdx++) {
            short[] roiImg = rois[roiIdx].getPixelsXYZ();
            for (int pixelIdx = 0; pixelIdx < (numROIX * numROIY); pixelIdx++) {
                if (roiImg[pixelIdx] > 0) {
                    roiFieldData[pixelIdx] = regionCounter;
                }
            }
            regionCounter++;
        }
        // create field data
        int NumTimePoints = 1;
        // 8 rois integrated into 1 image
        int NumChannels = 1;
        short[][][] pixData = new short[NumTimePoints][NumChannels][];
        pixData[0][0] = roiFieldData;
        // get extental data id
        VCImage vcImage = new VCImageUncompressed(null, new byte[isize.getXYZ()], extent, isize.getX(), isize.getY(), isize.getZ());
        RegionImage regionImage = new RegionImage(vcImage, 0, null, null, RegionImage.NO_SMOOTHING);
        CartesianMesh simpleCartesianMesh = CartesianMesh.createSimpleCartesianMesh(origin, extent, isize, regionImage);
        ExternalDataIdentifier newROIExtDataID = createNewExternalDataInfo(localWorkspace, ROI_SUMDATA_NAME).getExternalDataIdentifier();
        try {
            FieldDataFileOperationSpec fdos = new FieldDataFileOperationSpec();
            fdos.opType = FieldDataFileOperationSpec.FDOS_ADD;
            fdos.cartesianMesh = simpleCartesianMesh;
            fdos.shortSpecData = pixData;
            fdos.specEDI = newROIExtDataID;
            fdos.varNames = new String[] { "roiSumDataVar" };
            fdos.owner = LocalWorkspace.getDefaultOwner();
            fdos.times = new double[] { 0.0 };
            fdos.variableTypes = new VariableType[] { VariableType.VOLUME };
            fdos.origin = origin;
            fdos.extent = extent;
            fdos.isize = isize;
            localWorkspace.getDataSetControllerImpl().fieldDataFileOperation(fdos);
        } catch (Exception e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
        return new ROIDataGenerator(ROI_EXTDATA_NAME, /*name*/
        new int[] { 0 }, /* volumePoints*/
        new int[0], /* membranePoints*/
        regionCounter, /*numRegions*/
        0, /*zSlice*/
        newROIExtDataID.getKey(), /* fieldDataKey, sample image*/
        new FieldFunctionArguments(ROI_SUMDATA_NAME, ROI_SUMDATA_VARNAME, new Expression(0), VariableType.VOLUME), /*FieldFunctionArguments, sample image*/
        false);
    }
    return null;
}
Also used : Origin(org.vcell.util.Origin) Extent(org.vcell.util.Extent) FieldFunctionArguments(cbit.vcell.field.FieldFunctionArguments) ISize(org.vcell.util.ISize) FieldDataFileOperationSpec(cbit.vcell.field.io.FieldDataFileOperationSpec) VCImage(cbit.image.VCImage) VCImageUncompressed(cbit.image.VCImageUncompressed) ImageException(cbit.image.ImageException) ObjectNotFoundException(org.vcell.util.ObjectNotFoundException) IOException(java.io.IOException) UserCancelException(org.vcell.util.UserCancelException) ROIDataGenerator(org.vcell.vmicro.workflow.data.ROIDataGenerator) CartesianMesh(cbit.vcell.solvers.CartesianMesh) Expression(cbit.vcell.parser.Expression) RegionImage(cbit.vcell.geometry.RegionImage) ExternalDataIdentifier(org.vcell.util.document.ExternalDataIdentifier)

Aggregations

Origin (org.vcell.util.Origin)64 Extent (org.vcell.util.Extent)58 ISize (org.vcell.util.ISize)45 CartesianMesh (cbit.vcell.solvers.CartesianMesh)18 VCImageUncompressed (cbit.image.VCImageUncompressed)17 FieldDataFileOperationSpec (cbit.vcell.field.io.FieldDataFileOperationSpec)17 ImageException (cbit.image.ImageException)14 RegionImage (cbit.vcell.geometry.RegionImage)14 VCImage (cbit.image.VCImage)13 ArrayList (java.util.ArrayList)13 Expression (cbit.vcell.parser.Expression)12 IOException (java.io.IOException)12 SubVolume (cbit.vcell.geometry.SubVolume)11 File (java.io.File)11 UShortImage (cbit.vcell.VirtualMicroscopy.UShortImage)10 Geometry (cbit.vcell.geometry.Geometry)9 ExternalDataIdentifier (org.vcell.util.document.ExternalDataIdentifier)8 ImageDataset (cbit.vcell.VirtualMicroscopy.ImageDataset)7 BioModel (cbit.vcell.biomodel.BioModel)7 AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)7