Examples with RandomAccess net.imglib2.RandomAccess used on opensource projects

Example 1 with RandomAccess

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

the class DeconstructGeometryCommand method run.

@Override
public void run() {
    // Crop to get a z-stack over time (remove channel dimension)
    long maxX = fluorData.max(fluorData.dimensionIndex(Axes.X));
    long maxY = fluorData.max(fluorData.dimensionIndex(Axes.Y));
    long maxZ = fluorData.max(fluorData.dimensionIndex(Axes.Z));
    long maxTime = fluorData.max(fluorData.dimensionIndex(Axes.TIME));
    Img fluorImg = fluorData.getImgPlus().getImg();
    FinalInterval intervals = Intervals.createMinMax(0, 0, 0, 0, 0, maxX, maxY, maxZ, 0, maxTime);
    RandomAccessibleInterval fluorImgCropped = ops.transform().crop(fluorImg, intervals, true);
    // Calculate scale factors
    double[] scaleFactors = { 1, 1, 1, 1 };
    for (int i = 0; i < geomData.numDimensions(); i++) {
        scaleFactors[i] = geomData.dimension(i) / (double) fluorImgCropped.dimension(i);
    }
    // Scale the fluorescence dataset to match the geometry
    NLinearInterpolatorFactory interpolatorFactory = new NLinearInterpolatorFactory();
    RandomAccessibleInterval fluorScaled = ops.transform().scale(fluorImgCropped, scaleFactors, interpolatorFactory);
    // Crop out the first slice of each z-stack in time series
    intervals = Intervals.createMinMax(0, 0, 0, 0, fluorScaled.dimension(0) - 1, fluorScaled.dimension(1) - 1, 0, fluorScaled.dimension(3) - 1);
    IntervalView fluorXYT = (IntervalView) ops.transform().crop(fluorScaled, intervals, true);
    // Create a blank image of the same X-Y-Time dimensions
    long[] dimensions = { fluorXYT.dimension(0), fluorXYT.dimension(1), fluorXYT.dimension(2) };
    Img<DoubleType> result = ops.create().img(dimensions);
    // Calculate constant d in TIRF exponential decay function
    theta = theta * 2 * Math.PI / 360;
    double n1 = 1.52;
    double n2 = 1.38;
    double d = lambda * Math.pow((Math.pow(n1, 2) * Math.pow(Math.sin(theta), 2) - Math.pow(n2, 2)), -0.5) / (4 * Math.PI);
    // Iterate through each time point, using 3D geometry to generate 2D intensities
    Cursor<DoubleType> cursor = fluorXYT.localizingCursor();
    RandomAccess fluorRA = fluorScaled.randomAccess();
    RandomAccess<RealType<?>> geomRA = geomData.randomAccess();
    RandomAccess<DoubleType> resultRA = result.randomAccess();
    maxZ = geomData.dimension(2) - 1;
    while (cursor.hasNext()) {
        cursor.fwd();
        int[] positionXYZ = { cursor.getIntPosition(0), cursor.getIntPosition(1), (int) maxZ - 1 };
        int[] positionXYZT = { cursor.getIntPosition(0), cursor.getIntPosition(1), (int) maxZ - 1, cursor.getIntPosition(2) };
        resultRA.setPosition(cursor);
        geomRA.setPosition(positionXYZ);
        double sum = 0.0;
        while (positionXYZ[2] >= 0 && geomRA.get().getRealDouble() != 0.0) {
            fluorRA.setPosition(positionXYZT);
            geomRA.setPosition(positionXYZ);
            sum += geomRA.get().getRealDouble() * Math.exp(-zSpacing * positionXYZ[2] / d);
            positionXYZ[2]--;
        }
        resultRA.get().set(sum);
    }
    System.out.println("done");
    displayService.createDisplay(result);
}

Example 2 with RandomAccess

use of net.imglib2.RandomAccess 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");
}

Example 3 with RandomAccess

use of net.imglib2.RandomAccess 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);
}

Example 4 with RandomAccess

use of net.imglib2.RandomAccess 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;
}