use of cbit.vcell.geometry.ImageSubVolume in project vcell by virtualcell.
the class RayCaster method resampleGeometry.
public static Geometry resampleGeometry(GeometryThumbnailImageFactory geometryThumbnailImageFactory, Geometry origGeometry, ISize sampleSize) throws ImageException, PropertyVetoException, GeometryException, ExpressionException {
if (origGeometry.getDimension() < 3) {
throw new GeometryException("Presently, the Raycaster resampling works only for 3d geometries.");
}
GeometrySpec origGeometrySpec = origGeometry.getGeometrySpec();
VCImage origSubvolumeImage = origGeometrySpec.getSampledImage().getCurrentValue();
if (origSubvolumeImage == null) {
throw new GeometryException("original geometry does not have a sampled image");
}
VCImage resampledSubvolumeImage = RayCaster.sampleGeometry(origGeometry, sampleSize, false);
//
// Check if resampling failed:
// if not the same number of pixelClasses (between original geometry and resampled)
// if the subvolume handles are the same
//
boolean bSameSubvolumes = true;
if (origSubvolumeImage.getNumPixelClasses() != resampledSubvolumeImage.getNumPixelClasses()) {
bSameSubvolumes = false;
}
//
for (VCPixelClass origPixelClass : origSubvolumeImage.getPixelClasses()) {
VCPixelClass resampledPixelClass = resampledSubvolumeImage.getPixelClassFromPixelValue(origPixelClass.getPixel());
if (resampledPixelClass == null) {
bSameSubvolumes = false;
break;
}
}
//
if (!bSameSubvolumes) {
StringBuffer message = new StringBuffer();
message.append("\n\nexisting geometry:\n");
for (SubVolume oldSubvolume : origGeometrySpec.getSubVolumes()) {
long count = origSubvolumeImage.countPixelsByValue((byte) oldSubvolume.getHandle());
message.append("subvolume('" + oldSubvolume.getName() + "',handle=" + oldSubvolume.getHandle() + ",numPixels=" + count + " of " + origSubvolumeImage.getNumXYZ() + "\n");
}
message.append("\n\nnew resampled handle VCImage:\n");
for (VCPixelClass newPixelClass : resampledSubvolumeImage.getPixelClasses()) {
long count = resampledSubvolumeImage.countPixelsByValue((byte) newPixelClass.getPixel());
message.append("pixelClass('" + newPixelClass.getPixelClassName() + "',pixelValue=" + newPixelClass.getPixel() + ",numPixels=" + count + " of " + resampledSubvolumeImage.getNumXYZ() + ")\n");
}
throw new GeometryException("original Geometry had " + origSubvolumeImage.getNumPixelClasses() + " subvolumes, resampled Geometry found " + resampledSubvolumeImage.getNumPixelClasses() + " subvolumes " + message.toString());
}
//
// Create new VCImage that will form the basis for a new image-based geometry.
//
VCImage newVCImage = null;
if (origGeometrySpec.getImage() != null) {
//
// was an image-based geometry - try to make new VCImage similar to the original (same pixelClass names and pixel values).
// the goal is to make identical geometries if the sample size is same as the original image size.
//
// create a new VCImage with same image pixel values (not subvolume handles) and pixel class names as original image.
//
byte[] newVCImagePixels = new byte[sampleSize.getXYZ()];
byte[] resampledSubvolumePixels = resampledSubvolumeImage.getPixels();
for (int i = 0; i < sampleSize.getXYZ(); i++) {
int subvolumeHandle = resampledSubvolumePixels[i];
ImageSubVolume imageSubvolume = (ImageSubVolume) origGeometrySpec.getSubVolume(subvolumeHandle);
newVCImagePixels[i] = (byte) imageSubvolume.getPixelValue();
}
newVCImage = new VCImageUncompressed(null, newVCImagePixels, origGeometry.getExtent(), sampleSize.getX(), sampleSize.getY(), sampleSize.getZ());
newVCImage.setName(origGeometrySpec.getImage().getName());
ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
for (VCPixelClass origPixelClass : origGeometrySpec.getImage().getPixelClasses()) {
SubVolume origSubvolume = origGeometrySpec.getImageSubVolumeFromPixelValue(origPixelClass.getPixel());
newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origPixelClass.getPixel()));
}
newVCImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
} else {
//
// was an analytic geometry - create a new image-based geometry
//
// create a new VCImage with image pixel values and pixelClass names equal to corresponding subvolumes from original geometry
// make new subvolume names equal to old subvolume names.
//
byte[] newVCImageSubvolumePixels = resampledSubvolumeImage.getPixels().clone();
newVCImage = new VCImageUncompressed(null, newVCImageSubvolumePixels, origGeometry.getExtent(), sampleSize.getX(), sampleSize.getY(), sampleSize.getZ());
ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
for (SubVolume origSubvolume : origGeometrySpec.getSubVolumes()) {
// note: newVCImage already has subvolume handle pixels.
newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origSubvolume.getHandle()));
}
newVCImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
}
//
// construct the new geometry with the sampled VCImage.
//
Geometry newGeometry = new Geometry(origGeometry.getName(), newVCImage);
newGeometry.getGeometrySpec().setExtent(origGeometry.getExtent());
newGeometry.getGeometrySpec().setOrigin(origGeometry.getOrigin());
newGeometry.setDescription(origGeometry.getDescription());
newGeometry.getGeometrySurfaceDescription().setFilterCutoffFrequency(origGeometry.getGeometrySurfaceDescription().getFilterCutoffFrequency());
newGeometry.precomputeAll(geometryThumbnailImageFactory, true, true);
return newGeometry;
}
use of cbit.vcell.geometry.ImageSubVolume in project vcell by virtualcell.
the class GeometrySubVolumePanel method refreshButtons.
/**
* This method was created in VisualAge.
*/
private void refreshButtons() {
boolean bHasGeometry = getGeometry() != null;
boolean bHasGeomSpec = (bHasGeometry ? getGeometry().getGeometrySpec() != null : false);
boolean bImageBased = (bHasGeomSpec ? getGeometry().getGeometrySpec().getImage() != null : false);
addShapeButton.setEnabled(!bImageBased);
SubVolume selectedSubVolume = getSelectedSubVolume();
if (!bHasGeomSpec || getGeometry().getDimension() == 0) {
getFrontButton().setEnabled(false);
getBackButton().setEnabled(false);
getDeleteButton().setEnabled(false);
addShapeButton.setEnabled(false);
} else if (selectedSubVolume == null) {
getFrontButton().setEnabled(false);
getBackButton().setEnabled(false);
getDeleteButton().setEnabled(false);
} else {
GeometrySpec geometrySpec = getGeometry().getGeometrySpec();
int numAnalyticOrCSGSubVolumes = geometrySpec.getNumAnalyticOrCSGSubVolumes();
if (numAnalyticOrCSGSubVolumes > 1) {
getFrontButton().setEnabled(geometrySpec.getSubVolumeIndex(selectedSubVolume) > 0);
getBackButton().setEnabled(geometrySpec.getSubVolumeIndex(selectedSubVolume) < (numAnalyticOrCSGSubVolumes - 1));
} else {
getFrontButton().setEnabled(false);
getBackButton().setEnabled(false);
}
if (selectedSubVolume instanceof ImageSubVolume || (!bImageBased && (geometrySpec.getNumSubVolumes() <= 1))) {
getDeleteButton().setEnabled(false);
} else {
getDeleteButton().setEnabled(true);
}
}
}
use of cbit.vcell.geometry.ImageSubVolume 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.ImageSubVolume 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
// }
}
use of cbit.vcell.geometry.ImageSubVolume in project vcell by virtualcell.
the class SubVolumeTable method getImageSubVolume.
/**
* This method was created in VisualAge.
* @return Model
* @param rset ResultSet
* @param log SessionLog
*/
public ImageSubVolume getImageSubVolume(KeyValue key, ResultSet rset, cbit.image.VCPixelClass pixelClass) throws SQLException, DataAccessException {
// KeyValue key = new KeyValue(rset.getBigDecimal(id.toString(),0));
String svName = rset.getString(name.toString());
int handleValue = rset.getInt(handle.toString());
ImageSubVolume imageSubVolume = new ImageSubVolume(key, pixelClass, handleValue);
try {
imageSubVolume.setName(svName);
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new DataAccessException(e.getMessage());
}
return imageSubVolume;
}
Aggregations