use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class XmlReader method getMembraneSubDomain.
/**
* This method returns a MembraneSubDomain object from a XML Element.
* Creation date: (5/18/2001 4:23:30 PM)
* @return cbit.vcell.math.MembraneSubDomain
* @param param org.jdom.Element
* @exception cbit.vcell.xml.XmlParseException The exception description.
*/
@SuppressWarnings("unchecked")
private MembraneSubDomain getMembraneSubDomain(Element param, MathDescription mathDesc) throws XmlParseException {
// no need to do anything with the 'Name' attribute : constructor of MembraneSubDomain creates name from inside/outside compartmentSubDomains.
// String msdName = unMangle( param.getAttributeValue(XMLTags.NameAttrTag) );
// if ( msdName != null) {
// }
// get compartmentSubDomain references
// inside
String name = unMangle(param.getAttributeValue(XMLTags.InsideCompartmentTag));
CompartmentSubDomain insideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
if (insideRef == null) {
throw new XmlParseException("The reference to the inside CompartmentSubDomain " + name + ", could not be resolved!");
}
// outside
name = unMangle(param.getAttributeValue(XMLTags.OutsideCompartmentTag));
CompartmentSubDomain outsideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
if (outsideRef == null) {
throw new XmlParseException("The reference to the outside CompartmentSubDomain " + name + ", could not be resolved!");
}
// *** create new Membrane SubDomain ***
SubVolume insideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(insideRef.getName());
SubVolume outsideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(outsideRef.getName());
SurfaceClass surfaceClass = mathDesc.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubVolume, outsideSubVolume);
MembraneSubDomain subDomain = new MembraneSubDomain(insideRef, outsideRef, surfaceClass.getName());
transcribeComments(param, subDomain);
// Process BoundaryConditions
Iterator<Element> iterator = param.getChildren(XMLTags.BoundaryTypeTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
// create BoundaryConditionType
String temp = tempelement.getAttributeValue(XMLTags.BoundaryTypeAttrTag);
BoundaryConditionType bType = new BoundaryConditionType(temp);
// Process Xm
if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXm)) {
subDomain.setBoundaryConditionXm(bType);
} else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXp)) {
// Process Xp
subDomain.setBoundaryConditionXp(bType);
} else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYm)) {
// Process Ym
subDomain.setBoundaryConditionYm(bType);
} else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYp)) {
// Process Yp
subDomain.setBoundaryConditionYp(bType);
} else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZm)) {
// Process Zm
subDomain.setBoundaryConditionZm(bType);
} else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZp)) {
// Process Zp
subDomain.setBoundaryConditionZp(bType);
} else {
// If not indentified throw an exception!!
throw new XmlParseException("Unknown BoundaryConditionType: " + tempelement.getAttributeValue(XMLTags.BoundaryAttrTag));
}
}
// Add OdeEquations
iterator = param.getChildren(XMLTags.OdeEquationTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempElement = (Element) iterator.next();
OdeEquation odeEquation = getOdeEquation(tempElement, mathDesc);
try {
subDomain.addEquation(odeEquation);
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding an OdeEquation to a MembraneSubDomain!", e);
}
}
// process PdeEquations
iterator = param.getChildren(XMLTags.PdeEquationTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempElement = (Element) iterator.next();
try {
subDomain.addEquation(getPdeEquation(tempElement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding an PdeEquation to the MembraneSubDomain " + name, e);
}
}
// Add JumpConditions
iterator = param.getChildren(XMLTags.JumpConditionTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempElement = (Element) iterator.next();
try {
subDomain.addJumpCondition(getJumpCondition(tempElement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding a JumpCondition to a MembraneSubDomain!", e);
}
}
// Add the FastSystem (if any)
Element tempElement = param.getChild(XMLTags.FastSystemTag, vcNamespace);
if (tempElement != null) {
subDomain.setFastSystem(getFastSystem(tempElement, mathDesc));
}
// add MembraneRegionEquation
iterator = param.getChildren(XMLTags.MembraneRegionEquationTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempElement = (Element) iterator.next();
try {
subDomain.addEquation(getMembraneRegionEquation(tempElement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding a MembraneRegionEquation to a MEmbraneSubDomain!", e);
}
}
iterator = param.getChildren(XMLTags.ParticleJumpProcessTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
try {
subDomain.addParticleJumpProcess(getParticleJumpProcess(tempelement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
}
}
iterator = param.getChildren(XMLTags.ParticlePropertiesTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
try {
subDomain.addParticleProperties(getParticleProperties(tempelement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
}
}
// process ComputeNormal "equations"
iterator = param.getChildren(XMLTags.ComputeNormalTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
try {
subDomain.addEquation(getComputeNormal(tempelement, mathDesc));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("A MathException was fired when adding an ComputeNormal 'equation' to the MembraneSubDomain " + name, e);
}
}
Element velElem = param.getChild(XMLTags.VelocityTag, vcNamespace);
setMembraneSubdomainVelocity(velElem, XMLTags.XAttrTag, subDomain::setVelocityX);
setMembraneSubdomainVelocity(velElem, XMLTags.YAttrTag, subDomain::setVelocityY);
return subDomain;
}
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class XmlReader method getSurfaceClass.
private SurfaceClass getSurfaceClass(Element param, Geometry geom) throws XmlParseException {
Set<SubVolume> surfaceClassSubVolumeSet = new HashSet<SubVolume>();
KeyValue surfaceClassKey = null;
String surfaceClassName = null;
surfaceClassName = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
String surfaceClassKeyStr = param.getAttributeValue(XMLTags.KeyValueAttrTag);
surfaceClassKey = (surfaceClassKeyStr == null ? null : new KeyValue(surfaceClassKeyStr));
String subVol1Ref = param.getAttributeValue(XMLTags.SubVolume1RefAttrTag);
String subVol2Ref = param.getAttributeValue(XMLTags.SubVolume2RefAttrTag);
if (subVol1Ref != null) {
SubVolume subVolume = geom.getGeometrySpec().getSubVolume(subVol1Ref);
if (subVolume == null) {
throw new XmlParseException("SurfaceClass missing subvolume '" + subVol1Ref + "'");
}
surfaceClassSubVolumeSet.add(subVolume);
}
if (subVol2Ref != null) {
SubVolume subVolume = geom.getGeometrySpec().getSubVolume(subVol2Ref);
if (subVolume == null) {
throw new XmlParseException("SurfaceClass missing subvolume '" + subVol2Ref + "'");
}
surfaceClassSubVolumeSet.add(subVolume);
}
return new SurfaceClass(surfaceClassSubVolumeSet, surfaceClassKey, surfaceClassName);
}
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class OutputFunctionsPanel method addFunction.
private void addFunction() throws Exception {
String funcName = getFunctionNameTextField().getText();
Expression funcExp = null;
funcExp = new Expression(getFunctionExpressionTextField().getText());
Domain domain = null;
VariableType newFunctionVariableType = null;
boolean bSpatial = simulationWorkspace.getSimulationOwner().getGeometry().getDimension() > 0;
if (bSpatial) {
Object selectedItem = getSubdomainComboBox().getSelectedItem();
if (selectedItem instanceof GeometryClass) {
GeometryClass geoClass = (GeometryClass) selectedItem;
domain = new Domain(geoClass);
if (selectedItem instanceof SubVolume) {
newFunctionVariableType = VariableType.VOLUME;
} else if (selectedItem instanceof SurfaceClass) {
newFunctionVariableType = VariableType.MEMBRANE;
}
} else if (selectedItem instanceof VariableType) {
newFunctionVariableType = (VariableType) selectedItem;
} else {
newFunctionVariableType = VariableType.UNKNOWN;
}
} else {
newFunctionVariableType = VariableType.NONSPATIAL;
}
AnnotatedFunction tempFunction = new AnnotatedFunction(funcName, funcExp, domain, null, newFunctionVariableType, FunctionCategory.OUTPUTFUNCTION);
VariableType vt = outputFunctionContext.computeFunctionTypeWRTExpression(tempFunction, funcExp);
FunctionCategory category = FunctionCategory.OUTPUTFUNCTION;
if (vt.equals(VariableType.POSTPROCESSING)) {
category = FunctionCategory.POSTPROCESSFUNCTION;
}
AnnotatedFunction newFunction = new AnnotatedFunction(funcName, funcExp, domain, null, vt, category);
outputFunctionContext.addOutputFunction(newFunction);
setSelectedObjects(new Object[] { newFunction });
enableDeleteFnButton();
}
use of cbit.vcell.geometry.SurfaceClass in project vcell by virtualcell.
the class FRAPStudy method createNewSimBioModel.
public static BioModel createNewSimBioModel(FRAPStudy sourceFrapStudy, Parameter[] params, TimeStep tStep, KeyValue simKey, User owner, int startingIndexForRecovery) throws Exception {
if (owner == null) {
throw new Exception("Owner is not defined");
}
ROI cellROI_2D = sourceFrapStudy.getFrapData().getRoi(FRAPData.VFRAP_ROI_ENUM.ROI_CELL.name());
double df = params[FRAPModel.INDEX_PRIMARY_DIFF_RATE].getInitialGuess();
double ff = params[FRAPModel.INDEX_PRIMARY_FRACTION].getInitialGuess();
double bwmRate = params[FRAPModel.INDEX_BLEACH_MONITOR_RATE].getInitialGuess();
double dc = 0;
double fc = 0;
double bs = 0;
double onRate = 0;
double offRate = 0;
if (params.length == FRAPModel.NUM_MODEL_PARAMETERS_TWO_DIFF) {
dc = params[FRAPModel.INDEX_SECONDARY_DIFF_RATE].getInitialGuess();
fc = params[FRAPModel.INDEX_SECONDARY_FRACTION].getInitialGuess();
} else if (params.length == FRAPModel.NUM_MODEL_PARAMETERS_BINDING) {
dc = params[FRAPModel.INDEX_SECONDARY_DIFF_RATE].getInitialGuess();
fc = params[FRAPModel.INDEX_SECONDARY_FRACTION].getInitialGuess();
bs = params[FRAPModel.INDEX_BINDING_SITE_CONCENTRATION].getInitialGuess();
onRate = params[FRAPModel.INDEX_ON_RATE].getInitialGuess();
offRate = params[FRAPModel.INDEX_OFF_RATE].getInitialGuess();
}
// immobile fraction
double fimm = 1 - ff - fc;
if (fimm < FRAPOptimizationUtils.epsilon && fimm > (0 - FRAPOptimizationUtils.epsilon)) {
fimm = 0;
}
if (fimm < (1 + FRAPOptimizationUtils.epsilon) && fimm > (1 - FRAPOptimizationUtils.epsilon)) {
fimm = 1;
}
Extent extent = sourceFrapStudy.getFrapData().getImageDataset().getExtent();
double[] timeStamps = sourceFrapStudy.getFrapData().getImageDataset().getImageTimeStamps();
TimeBounds timeBounds = new TimeBounds(0.0, timeStamps[timeStamps.length - 1] - timeStamps[startingIndexForRecovery]);
double timeStepVal = timeStamps[startingIndexForRecovery + 1] - timeStamps[startingIndexForRecovery];
int numX = cellROI_2D.getRoiImages()[0].getNumX();
int numY = cellROI_2D.getRoiImages()[0].getNumY();
int numZ = cellROI_2D.getRoiImages().length;
short[] shortPixels = cellROI_2D.getRoiImages()[0].getPixels();
byte[] bytePixels = new byte[numX * numY * numZ];
final byte EXTRACELLULAR_PIXVAL = 0;
final byte CYTOSOL_PIXVAL = 1;
for (int i = 0; i < bytePixels.length; i++) {
if (shortPixels[i] != 0) {
bytePixels[i] = CYTOSOL_PIXVAL;
}
}
VCImage maskImage;
try {
maskImage = new VCImageUncompressed(null, bytePixels, extent, numX, numY, numZ);
} catch (ImageException e) {
e.printStackTrace();
throw new RuntimeException("failed to create mask image for geometry");
}
Geometry geometry = new Geometry("geometry", maskImage);
if (geometry.getGeometrySpec().getNumSubVolumes() != 2) {
throw new Exception("Cell ROI has no ExtraCellular.");
}
int subVolume0PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(0)).getPixelValue();
geometry.getGeometrySpec().getSubVolume(0).setName((subVolume0PixVal == EXTRACELLULAR_PIXVAL ? EXTRACELLULAR_NAME : CYTOSOL_NAME));
int subVolume1PixVal = ((ImageSubVolume) geometry.getGeometrySpec().getSubVolume(1)).getPixelValue();
geometry.getGeometrySpec().getSubVolume(1).setName((subVolume1PixVal == CYTOSOL_PIXVAL ? CYTOSOL_NAME : EXTRACELLULAR_NAME));
geometry.getGeometrySurfaceDescription().updateAll();
BioModel bioModel = new BioModel(null);
bioModel.setName("unnamed");
Model model = new Model("model");
bioModel.setModel(model);
model.addFeature(EXTRACELLULAR_NAME);
Feature extracellular = (Feature) model.getStructure(EXTRACELLULAR_NAME);
model.addFeature(CYTOSOL_NAME);
Feature cytosol = (Feature) model.getStructure(CYTOSOL_NAME);
// Membrane mem = model.addMembrane(EXTRACELLULAR_CYTOSOL_MEM_NAME);
// model.getStructureTopology().setInsideFeature(mem, cytosol);
// model.getStructureTopology().setOutsideFeature(mem, extracellular);
String roiDataName = FRAPStudy.ROI_EXTDATA_NAME;
final int SPECIES_COUNT = 4;
final int FREE_SPECIES_INDEX = 0;
final int BS_SPECIES_INDEX = 1;
final int COMPLEX_SPECIES_INDEX = 2;
final int IMMOBILE_SPECIES_INDEX = 3;
Expression[] diffusionConstants = null;
Species[] species = null;
SpeciesContext[] speciesContexts = null;
Expression[] initialConditions = null;
diffusionConstants = new Expression[SPECIES_COUNT];
species = new Species[SPECIES_COUNT];
speciesContexts = new SpeciesContext[SPECIES_COUNT];
initialConditions = new Expression[SPECIES_COUNT];
// total initial condition
FieldFunctionArguments postBleach_first = new FieldFunctionArguments(roiDataName, "postbleach_first", new Expression(0), VariableType.VOLUME);
FieldFunctionArguments prebleach_avg = new FieldFunctionArguments(roiDataName, "prebleach_avg", new Expression(0), VariableType.VOLUME);
Expression expPostBleach_first = new Expression(postBleach_first.infix());
Expression expPreBleach_avg = new Expression(prebleach_avg.infix());
Expression totalIniCondition = Expression.div(expPostBleach_first, expPreBleach_avg);
// Free Species
diffusionConstants[FREE_SPECIES_INDEX] = new Expression(df);
species[FREE_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_MOBILE, "Mobile bleachable species");
speciesContexts[FREE_SPECIES_INDEX] = new SpeciesContext(null, species[FREE_SPECIES_INDEX].getCommonName(), species[FREE_SPECIES_INDEX], cytosol);
initialConditions[FREE_SPECIES_INDEX] = Expression.mult(new Expression(ff), totalIniCondition);
// Immobile Species (No diffusion)
// Set very small diffusion rate on immobile to force evaluation as state variable (instead of FieldData function)
// If left as a function errors occur because functions involving FieldData require a database connection
final String IMMOBILE_DIFFUSION_KLUDGE = "1e-14";
diffusionConstants[IMMOBILE_SPECIES_INDEX] = new Expression(IMMOBILE_DIFFUSION_KLUDGE);
species[IMMOBILE_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_IMMOBILE, "Immobile bleachable species");
speciesContexts[IMMOBILE_SPECIES_INDEX] = new SpeciesContext(null, species[IMMOBILE_SPECIES_INDEX].getCommonName(), species[IMMOBILE_SPECIES_INDEX], cytosol);
initialConditions[IMMOBILE_SPECIES_INDEX] = Expression.mult(new Expression(fimm), totalIniCondition);
// BS Species
diffusionConstants[BS_SPECIES_INDEX] = new Expression(IMMOBILE_DIFFUSION_KLUDGE);
species[BS_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_BINDING_SITE, "Binding Site species");
speciesContexts[BS_SPECIES_INDEX] = new SpeciesContext(null, species[BS_SPECIES_INDEX].getCommonName(), species[BS_SPECIES_INDEX], cytosol);
initialConditions[BS_SPECIES_INDEX] = Expression.mult(new Expression(bs), totalIniCondition);
// Complex species
diffusionConstants[COMPLEX_SPECIES_INDEX] = new Expression(dc);
species[COMPLEX_SPECIES_INDEX] = new Species(FRAPStudy.SPECIES_NAME_PREFIX_SLOW_MOBILE, "Slower mobile bleachable species");
speciesContexts[COMPLEX_SPECIES_INDEX] = new SpeciesContext(null, species[COMPLEX_SPECIES_INDEX].getCommonName(), species[COMPLEX_SPECIES_INDEX], cytosol);
initialConditions[COMPLEX_SPECIES_INDEX] = Expression.mult(new Expression(fc), totalIniCondition);
// add reactions to species if there is bleachWhileMonitoring rate.
for (int i = 0; i < initialConditions.length; i++) {
model.addSpecies(species[i]);
model.addSpeciesContext(speciesContexts[i]);
// reaction with BMW rate, which should not be applied to binding site
if (!(species[i].getCommonName().equals(FRAPStudy.SPECIES_NAME_PREFIX_BINDING_SITE))) {
SimpleReaction simpleReaction = new SimpleReaction(model, cytosol, speciesContexts[i].getName() + "_bleach", true);
model.addReactionStep(simpleReaction);
simpleReaction.addReactant(speciesContexts[i], 1);
MassActionKinetics massActionKinetics = new MassActionKinetics(simpleReaction);
simpleReaction.setKinetics(massActionKinetics);
KineticsParameter kforward = massActionKinetics.getForwardRateParameter();
simpleReaction.getKinetics().setParameterValue(kforward, new Expression(new Double(bwmRate)));
}
}
// add the binding reaction: F + BS <-> C
SimpleReaction simpleReaction2 = new SimpleReaction(model, cytosol, "reac_binding", true);
model.addReactionStep(simpleReaction2);
simpleReaction2.addReactant(speciesContexts[FREE_SPECIES_INDEX], 1);
simpleReaction2.addReactant(speciesContexts[BS_SPECIES_INDEX], 1);
simpleReaction2.addProduct(speciesContexts[COMPLEX_SPECIES_INDEX], 1);
MassActionKinetics massActionKinetics = new MassActionKinetics(simpleReaction2);
simpleReaction2.setKinetics(massActionKinetics);
KineticsParameter kforward = massActionKinetics.getForwardRateParameter();
KineticsParameter kreverse = massActionKinetics.getReverseRateParameter();
simpleReaction2.getKinetics().setParameterValue(kforward, new Expression(new Double(onRate)));
simpleReaction2.getKinetics().setParameterValue(kreverse, new Expression(new Double(offRate)));
// create simulation context
SimulationContext simContext = new SimulationContext(bioModel.getModel(), geometry);
bioModel.addSimulationContext(simContext);
FeatureMapping cytosolFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(cytosol);
FeatureMapping extracellularFeatureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(extracellular);
// Membrane plasmaMembrane = model.getStructureTopology().getMembrane(cytosol, extracellular);
// MembraneMapping plasmaMembraneMapping = (MembraneMapping)simContext.getGeometryContext().getStructureMapping(plasmaMembrane);
SubVolume cytSubVolume = geometry.getGeometrySpec().getSubVolume(CYTOSOL_NAME);
SubVolume exSubVolume = geometry.getGeometrySpec().getSubVolume(EXTRACELLULAR_NAME);
SurfaceClass pmSurfaceClass = geometry.getGeometrySurfaceDescription().getSurfaceClass(exSubVolume, cytSubVolume);
cytosolFeatureMapping.setGeometryClass(cytSubVolume);
extracellularFeatureMapping.setGeometryClass(exSubVolume);
// plasmaMembraneMapping.setGeometryClass(pmSurfaceClass);
cytosolFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
extracellularFeatureMapping.getUnitSizeParameter().setExpression(new Expression(1.0));
for (int i = 0; i < speciesContexts.length; i++) {
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i]);
scs.getInitialConditionParameter().setExpression(initialConditions[i]);
scs.getDiffusionParameter().setExpression(diffusionConstants[i]);
}
MathMapping mathMapping = simContext.createNewMathMapping();
MathDescription mathDesc = mathMapping.getMathDescription();
// Add total fluorescence as function of mobile(optional: and slower mobile) and immobile fractions
mathDesc.addVariable(new Function(FRAPStudy.SPECIES_NAME_PREFIX_COMBINED, new Expression(species[FREE_SPECIES_INDEX].getCommonName() + "+" + species[COMPLEX_SPECIES_INDEX].getCommonName() + "+" + species[IMMOBILE_SPECIES_INDEX].getCommonName()), null));
simContext.setMathDescription(mathDesc);
SimulationVersion simVersion = new SimulationVersion(simKey, "sim1", owner, new GroupAccessNone(), new KeyValue("0"), new BigDecimal(0), new Date(), VersionFlag.Current, "", null);
Simulation newSimulation = new Simulation(simVersion, mathDesc);
simContext.addSimulation(newSimulation);
newSimulation.getSolverTaskDescription().setTimeBounds(timeBounds);
newSimulation.getMeshSpecification().setSamplingSize(cellROI_2D.getISize());
// newSimulation.getSolverTaskDescription().setTimeStep(timeStep); // Sundials doesn't need time step
newSimulation.getSolverTaskDescription().setSolverDescription(SolverDescription.SundialsPDE);
// use exp time step as output time spec
newSimulation.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec(timeStepVal));
return bioModel;
}
use of cbit.vcell.geometry.SurfaceClass 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
// }
}
Aggregations