Search in sources :

Example 6 with RandomAccessibleInterval

use of net.imglib2.RandomAccessibleInterval in project vcell by virtualcell.

the class ConstructTIRFGeometry method run.

@Override
public void run() {
    // Calculate constant d in TIRF exponential decay function
    // Angle of incidence in radians
    theta = theta * 2 * Math.PI / 360;
    // Refractive index of glass
    final double n1 = 1.52;
    // Refractive index of cytosol
    final double n2 = 1.38;
    final double d = lambda * Math.pow((Math.pow(n1, 2) * Math.pow(Math.sin(theta), 2) - Math.pow(n2, 2)), -0.5) / (4 * Math.PI);
    System.out.println("d: " + d);
    final double fluorPerMolecule = 250;
    // Get frame of interest to define geometry
    long maxX = data.dimension(0) - 1;
    long maxY = data.dimension(1) - 1;
    Interval interval = Intervals.createMinMax(0, 0, sliceIndex, maxX, maxY, sliceIndex);
    RandomAccessibleInterval<T> croppedRAI = ops.transform().crop(data, interval, true);
    // Subtract lowest pixel value
    IterableInterval<T> dataII = Views.iterable(croppedRAI);
    double min = ops.stats().min(dataII).getRealDouble();
    Cursor<T> dataCursor = dataII.cursor();
    while (dataCursor.hasNext()) {
        double val = dataCursor.next().getRealDouble();
        dataCursor.get().setReal(val - min);
    }
    // Perform Gaussian blur
    RandomAccessibleInterval<T> blurredRAI = ops.filter().gauss(croppedRAI, 2);
    IterableInterval<T> blurredII = Views.iterable(blurredRAI);
    // Segment slice by threshold and fill holes
    IterableInterval<BitType> thresholded = ops.threshold().huang(blurredII);
    Img<BitType> thresholdedImg = ops.convert().bit(thresholded);
    RandomAccessibleInterval<BitType> thresholdedRAI = ops.morphology().fillHoles(thresholdedImg);
    // Get the largest region
    RandomAccessibleInterval<LabelingType<ByteType>> labeling = ops.labeling().cca(thresholdedRAI, ConnectedComponents.StructuringElement.EIGHT_CONNECTED);
    LabelRegions<ByteType> labelRegions = new LabelRegions<>(labeling);
    Iterator<LabelRegion<ByteType>> iterator = labelRegions.iterator();
    LabelRegion<ByteType> maxRegion = iterator.next();
    while (iterator.hasNext()) {
        LabelRegion<ByteType> currRegion = iterator.next();
        if (currRegion.size() > maxRegion.size()) {
            maxRegion = currRegion;
        }
    }
    // Generate z index map
    double iMax = ops.stats().max(dataII).getRealDouble();
    Img<UnsignedShortType> dataImg = ops.convert().uint16(dataII);
    Img<UnsignedShortType> zMap = ops.convert().uint16(ops.create().img(dataII));
    LabelRegionCursor cursor = maxRegion.localizingCursor();
    RandomAccess<UnsignedShortType> zMapRA = zMap.randomAccess();
    RandomAccess<UnsignedShortType> dataRA = dataImg.randomAccess();
    while (cursor.hasNext()) {
        cursor.fwd();
        zMapRA.setPosition(cursor);
        dataRA.setPosition(cursor);
        double val = dataRA.get().getRealDouble();
        // Log of 0 is undefined
        if (val < 1) {
            val = 1;
        }
        int z = (int) Math.round(-d * Math.log(val / iMax) / zRes);
        zMapRA.get().set(z);
    }
    System.out.println("6");
    // Use map to construct 3D geometry
    // Add 5 slices of padding on top
    int maxZ = (int) ops.stats().max(zMap).getRealDouble() + 5;
    long[] resultDimensions = { maxX + 1, maxY + 1, maxZ };
    Img<BitType> result = new ArrayImgFactory<BitType>().create(resultDimensions, new BitType());
    RandomAccess<BitType> resultRA = result.randomAccess();
    System.out.println(maxZ);
    cursor.reset();
    while (cursor.hasNext()) {
        cursor.fwd();
        zMapRA.setPosition(cursor);
        int zIndex = zMapRA.get().get();
        int[] position = { cursor.getIntPosition(0), cursor.getIntPosition(1), zIndex };
        while (position[2] < maxZ) {
            resultRA.setPosition(position);
            resultRA.get().set(true);
            position[2]++;
        }
    }
    output = datasetService.create(result);
    CalibratedAxis[] axes = new DefaultLinearAxis[] { new DefaultLinearAxis(Axes.X), new DefaultLinearAxis(Axes.Y), new DefaultLinearAxis(Axes.Z) };
    output.setAxes(axes);
    System.out.println("Done constructing geometry");
}
Also used : ByteType(net.imglib2.type.numeric.integer.ByteType) DefaultLinearAxis(net.imagej.axis.DefaultLinearAxis) BitType(net.imglib2.type.logic.BitType) LabelRegionCursor(net.imglib2.roi.labeling.LabelRegionCursor) LabelingType(net.imglib2.roi.labeling.LabelingType) UnsignedShortType(net.imglib2.type.numeric.integer.UnsignedShortType) LabelRegion(net.imglib2.roi.labeling.LabelRegion) CalibratedAxis(net.imagej.axis.CalibratedAxis) LabelRegions(net.imglib2.roi.labeling.LabelRegions) RandomAccessibleInterval(net.imglib2.RandomAccessibleInterval) Interval(net.imglib2.Interval) IterableInterval(net.imglib2.IterableInterval)

Example 7 with RandomAccessibleInterval

use of net.imglib2.RandomAccessibleInterval in project vcell by virtualcell.

the class ImageStatsForPlotting method computeMean.

/**
 * Computes the mean of each XY slice along the 3rd dimension
 * TODO: Currently assumes only 3 dimensions, must handle time series of z stacks and multiple channels
 * @param data
 * @return Pair containing A) the 3rd dimension index, and B) the mean value of the XY slice
 */
private Pair<double[], double[]> computeMean(RandomAccessibleInterval<T> data, IterableInterval<BitType> mask) {
    double[] indices = new double[(int) data.dimension(2)];
    double[] means = new double[indices.length];
    for (int z = 0; z < indices.length; z++) {
        FinalInterval interval = Intervals.createMinMax(0, 0, z, data.dimension(0) - 1, data.dimension(1) - 1, z);
        double mean = 0.0;
        RandomAccessibleInterval<T> cropped = ops.transform().crop(data, interval);
        if (mask == null) {
            mean = ops.stats().mean(Views.iterable(cropped)).getRealDouble();
        } else {
            Cursor<BitType> maskCursor = mask.localizingCursor();
            RandomAccess<T> dataRA = cropped.randomAccess();
            RealSum sum = new RealSum();
            int size = 0;
            maskCursor.reset();
            while (maskCursor.hasNext()) {
                maskCursor.fwd();
                if (maskCursor.get().get()) {
                    dataRA.setPosition(maskCursor);
                    sum.add(dataRA.get().getRealDouble());
                    size++;
                }
            }
            mean = sum.getSum() / size;
        }
        indices[z] = z;
        means[z] = mean;
    }
    return new ValuePair<double[], double[]>(indices, means);
}
Also used : BitType(net.imglib2.type.logic.BitType) ValuePair(net.imglib2.util.ValuePair) FinalInterval(net.imglib2.FinalInterval) RealSum(net.imglib2.util.RealSum)

Example 8 with RandomAccessibleInterval

use of net.imglib2.RandomAccessibleInterval in project vcell by virtualcell.

the class LargestRegionSlice method run.

@Override
public void run() {
    int maxArea = 0;
    int zOfMaxArea = 0;
    for (int z = 0; z < data.dimension(2); z++) {
        FinalInterval interval = Intervals.createMinMax(0, 0, z, data.dimension(0) - 1, data.dimension(1) - 1, z);
        RandomAccessibleInterval<BitType> croppedRAI = ops.transform().crop(data, interval);
        IterableInterval<BitType> croppedII = Views.iterable(croppedRAI);
        Cursor<BitType> cursor = croppedII.cursor();
        int area = 0;
        while (cursor.hasNext()) {
            if (cursor.next().get()) {
                area++;
            }
        }
        if (area > maxArea) {
            maxArea = area;
            zOfMaxArea = z;
        }
    }
    FinalInterval interval = Intervals.createMinMax(0, 0, zOfMaxArea, data.dimension(0) - 1, data.dimension(1) - 1, zOfMaxArea);
    output = ops.transform().crop(data, interval);
}
Also used : BitType(net.imglib2.type.logic.BitType) FinalInterval(net.imglib2.FinalInterval)

Example 9 with RandomAccessibleInterval

use of net.imglib2.RandomAccessibleInterval in project vcell by virtualcell.

the class VCellService method runSimulation.

private static Task<List<Dataset>, SimulationState> runSimulation(final SimulationServiceImpl client, final VCellModel vCellModel, final SimulationSpec simSpec, final List<Species> outputSpecies, final boolean shouldCreateIndividualDatasets, final OpService opService, final DatasetService datasetService) throws IOException, XMLStreamException {
    final Task<List<Dataset>, SimulationState> task = new Task<List<Dataset>, SimulationState>() {

        @Override
        protected List<Dataset> doInBackground() throws Exception {
            setSubtask(SimulationState.notRun);
            final File sbmlSpatialFile = new File(vCellModel.getName() + ".xml");
            new SBMLWriter().write(vCellModel.getSbmlDocument(), sbmlSpatialFile);
            final SBMLModel model = new SBMLModel();
            model.setFilepath(sbmlSpatialFile.getAbsolutePath());
            final SimulationInfo simulationInfo = client.computeModel(model, simSpec);
            try {
                Thread.sleep(500);
            } catch (final InterruptedException e) {
                e.printStackTrace();
            }
            setSubtask(SimulationState.running);
            while (client.getStatus(simulationInfo).getSimState() == SimulationState.running) {
                System.out.println("waiting for simulation results");
                try {
                    Thread.sleep(500);
                } catch (final InterruptedException e) {
                    e.printStackTrace();
                }
            }
            if (client.getStatus(simulationInfo).getSimState() == SimulationState.failed) {
                setSubtask(SimulationState.failed);
                return null;
            }
            final List<Dataset> results = new ArrayList<>();
            final List<VariableInfo> vars = client.getVariableList(simulationInfo);
            final List<Double> times = client.getTimePoints(simulationInfo);
            for (final VariableInfo var : vars) {
                if (outputSpecies.stream().anyMatch(species -> species.getId().equals(var.getVariableVtuName()))) {
                    // Get data for first time point and determine dimensions
                    List<Double> data = client.getData(simulationInfo, var, 0);
                    final int[] dimensions = getDimensions(data, times);
                    final Img<DoubleType> img = opService.create().img(dimensions);
                    final RandomAccess<DoubleType> imgRA = img.randomAccess();
                    // Copy data to the ImgLib2 Img
                    for (int t = 0; t < times.size(); t++) {
                        data = client.getData(simulationInfo, var, t);
                        for (int d = 0; d < data.size(); d++) {
                            imgRA.setPosition(new int[] { d, t });
                            imgRA.get().set(data.get(d));
                        }
                    }
                    // Create ImageJ Dataset and add to results
                    final Dataset dataset = datasetService.create(img);
                    dataset.setName(var.getVariableVtuName());
                    results.add(dataset);
                }
            }
            // If desired, add all datasets with the same dimensions
            if (!shouldCreateIndividualDatasets && !results.isEmpty()) {
                // First, group datasets according to dimensions
                final List<List<Dataset>> datasetGroups = new ArrayList<>();
                final List<Dataset> initialGroup = new ArrayList<>();
                initialGroup.add(results.get(0));
                datasetGroups.add(initialGroup);
                for (int i = 1; i < results.size(); i++) {
                    final Dataset result = results.get(i);
                    for (final List<Dataset> datasetGroup : datasetGroups) {
                        final Dataset[] datasets = new Dataset[] { datasetGroup.get(0), result };
                        if (Datasets.areSameSize(datasets, 0, 1)) {
                            datasetGroup.add(result);
                        } else {
                            final List<Dataset> newGroup = new ArrayList<>();
                            newGroup.add(result);
                            datasetGroups.add(newGroup);
                        }
                    }
                }
                final List<Dataset> summedResults = new ArrayList<>();
                for (final List<Dataset> datasetGroup : datasetGroups) {
                    final Img<DoubleType> sum = opService.create().img(datasetGroup.get(0));
                    for (final Dataset dataset : datasetGroup) {
                        @SuppressWarnings("unchecked") final RandomAccessibleInterval<DoubleType> current = (Img<DoubleType>) dataset.getImgPlus().getImg();
                        opService.math().add(sum, sum, current);
                    }
                    final Dataset result = datasetService.create(sum);
                    result.setName(datasetGroup.stream().map(d -> d.getName()).collect(Collectors.joining("+")));
                    summedResults.add(result);
                }
                return summedResults;
            }
            setSubtask(SimulationState.done);
            return results;
        }
    };
    return task;
}
Also used : Task(org.vcell.imagej.common.gui.Task) SBMLModel(org.vcell.vcellij.api.SBMLModel) ArrayList(java.util.ArrayList) ArrayList(java.util.ArrayList) List(java.util.List) Img(net.imglib2.img.Img) Dataset(net.imagej.Dataset) VariableInfo(org.vcell.vcellij.api.VariableInfo) SBMLWriter(org.sbml.jsbml.SBMLWriter) DoubleType(net.imglib2.type.numeric.real.DoubleType) SimulationState(org.vcell.vcellij.api.SimulationState) File(java.io.File) SimulationInfo(org.vcell.vcellij.api.SimulationInfo)

Aggregations

Dataset (net.imagej.Dataset)4 FinalInterval (net.imglib2.FinalInterval)4 BitType (net.imglib2.type.logic.BitType)4 RandomAccessibleInterval (net.imglib2.RandomAccessibleInterval)3 Plot (ij.gui.Plot)2 Img (net.imglib2.img.Img)2 DoubleType (net.imglib2.type.numeric.real.DoubleType)2 Pair (net.imglib2.util.Pair)2 ColorPlot (org.vcell.imagej.common.gui.ColorPlot)2 File (java.io.File)1 ArrayList (java.util.ArrayList)1 List (java.util.List)1 CalibratedAxis (net.imagej.axis.CalibratedAxis)1 DefaultLinearAxis (net.imagej.axis.DefaultLinearAxis)1 Overlay (net.imagej.overlay.Overlay)1 Interval (net.imglib2.Interval)1 IterableInterval (net.imglib2.IterableInterval)1 RandomAccess (net.imglib2.RandomAccess)1 NLinearInterpolatorFactory (net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory)1 LabelRegion (net.imglib2.roi.labeling.LabelRegion)1