Examples with ImageStack ij.ImageStack used on opensource projects

Example 1 with ImageStack

use of ij.ImageStack in project GDSC-SMLM by aherbert.

the class PeakFit method run.

/**
	 * Locate the peaks in the configured image source. Results are saved to the configured output.
	 * <p>
	 * This must be called after initialisation with an image source. Note that each call to this method must be
	 * preceded with initialisation to prepare the image and output options.
	 */
public void run() {
    if (source == null)
        return;
    int totalFrames = source.getFrames();
    ImageStack stack = null;
    if (showProcessedFrames)
        stack = new ImageStack(bounds.width, bounds.height);
    // Do not crop the region from the source if the bounds match the source dimensions
    Rectangle cropBounds = (bounds.x == 0 && bounds.y == 0 && bounds.width == source.getWidth() && bounds.height == source.getHeight()) ? null : bounds;
    // Use the FitEngine to allow multi-threading.
    FitEngine engine = createFitEngine(getNumberOfThreads(totalFrames));
    final int step = Utils.getProgressInterval(totalFrames);
    runTime = System.nanoTime();
    boolean shutdown = false;
    int slice = 0;
    while (!shutdown) {
        // Noise can optionally be estimated from the entire frame
        float[] data = (ignoreBoundsForNoise) ? source.next() : source.next(cropBounds);
        if (data == null)
            break;
        if (++slice % step == 0) {
            if (Utils.showStatus("Slice: " + slice + " / " + totalFrames))
                IJ.showProgress(slice, totalFrames);
        }
        float noise = Float.NaN;
        if (ignoreBoundsForNoise) {
            noise = FitWorker.estimateNoise(data, source.getWidth(), source.getHeight(), config.getNoiseMethod());
            // Crop the data to the region
            data = ImageConverter.getData(data, source.getWidth(), source.getHeight(), bounds, null);
        }
        if (showProcessedFrames) {
            stack.addSlice(String.format("Frame %d - %d", source.getStartFrameNumber(), source.getEndFrameNumber()), data);
        }
        // Get the frame number from the source to allow for interlaced and aggregated data
        engine.run(createJob(source.getStartFrameNumber(), source.getEndFrameNumber(), data, bounds, noise));
        if (escapePressed())
            shutdown = true;
    }
    engine.end(shutdown);
    time = engine.getTime();
    runTime = System.nanoTime() - runTime;
    if (showProcessedFrames)
        Utils.display("Processed frames", stack);
    showResults();
    source.close();
}

Example 2 with ImageStack

use of ij.ImageStack in project GDSC-SMLM by aherbert.

the class DrawClusters method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog())
        return;
    // Load the results
    MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
    if (results == null || results.size() == 0) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    // Get the traces
    Trace[] traces = TraceManager.convert(results);
    if (traces == null || traces.length == 0) {
        IJ.error(TITLE, "No traces could be loaded");
        return;
    }
    // Filter traces to a min size
    int maxFrame = 0;
    int count = 0;
    final int myMaxSize = (maxSize < minSize) ? Integer.MAX_VALUE : maxSize;
    final boolean myDrawLines = (myMaxSize < 2) ? false : drawLines;
    for (int i = 0; i < traces.length; i++) {
        if (expandToSingles)
            traces[i].expandToSingles();
        if (traces[i].size() >= minSize && traces[i].size() <= myMaxSize) {
            traces[count++] = traces[i];
            traces[i].sort();
            if (maxFrame < traces[i].getTail().getFrame())
                maxFrame = traces[i].getTail().getFrame();
        }
    }
    if (count == 0) {
        IJ.error(TITLE, "No traces achieved the size limits");
        return;
    }
    String msg = String.format(TITLE + ": %d / %s (%s)", count, Utils.pleural(traces.length, "trace"), Utils.pleural(results.size(), "localisation"));
    IJ.showStatus(msg);
    //Utils.log(msg);
    Rectangle bounds = results.getBounds(true);
    ImagePlus imp = WindowManager.getImage(title);
    boolean isUseStackPosition = useStackPosition;
    if (imp == null) {
        // Create a default image using 100 pixels as the longest edge
        double maxD = (bounds.width > bounds.height) ? bounds.width : bounds.height;
        int w, h;
        if (maxD == 0) {
            // Note that imageSize can be zero (for auto sizing)
            w = h = (imageSize == 0) ? 20 : imageSize;
        } else {
            // Note that imageSize can be zero (for auto sizing)
            if (imageSize == 0) {
                w = bounds.width;
                h = bounds.height;
            } else {
                w = (int) (imageSize * bounds.width / maxD);
                h = (int) (imageSize * bounds.height / maxD);
            }
        }
        ByteProcessor bp = new ByteProcessor(w, h);
        if (isUseStackPosition) {
            ImageStack stack = new ImageStack(w, h, maxFrame);
            for (int i = 1; i <= maxFrame; i++) // Do not clone as the image is empty
            stack.setPixels(bp.getPixels(), i);
            imp = Utils.display(TITLE, stack);
        } else
            imp = Utils.display(TITLE, bp);
        // Enlarge
        ImageWindow iw = imp.getWindow();
        for (int i = 9; i-- > 0 && iw.getWidth() < 500 && iw.getHeight() < 500; ) {
            iw.getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
        }
    } else {
        // Check if the image has enough frames for all the traces
        if (maxFrame > imp.getNFrames())
            isUseStackPosition = false;
    }
    final float xScale = (float) (imp.getWidth() / bounds.getWidth());
    final float yScale = (float) (imp.getHeight() / bounds.getHeight());
    // Create ROIs and store data to sort them
    Roi[] rois = new Roi[count];
    int[][] frames = (isUseStackPosition) ? new int[count][] : null;
    int[] indices = Utils.newArray(count, 0, 1);
    double[] values = new double[count];
    for (int i = 0; i < count; i++) {
        Trace trace = traces[i];
        int nPoints = trace.size();
        float[] xPoints = new float[nPoints];
        float[] yPoints = new float[nPoints];
        int j = 0;
        if (isUseStackPosition)
            frames[i] = new int[nPoints];
        for (PeakResult result : trace.getPoints()) {
            xPoints[j] = (result.getXPosition() - bounds.x) * xScale;
            yPoints[j] = (result.getYPosition() - bounds.y) * yScale;
            if (isUseStackPosition)
                frames[i][j] = result.getFrame();
            j++;
        }
        Roi roi;
        if (myDrawLines) {
            roi = new PolygonRoi(xPoints, yPoints, nPoints, Roi.POLYLINE);
            if (splineFit)
                ((PolygonRoi) roi).fitSpline();
        } else {
            roi = new PointRoi(xPoints, yPoints, nPoints);
            ((PointRoi) roi).setShowLabels(false);
        }
        rois[i] = roi;
        switch(sort) {
            case 0:
            default:
                break;
            case // Sort by ID
            1:
                values[i] = traces[i].getId();
                break;
            case // Sort by time
            2:
                values[i] = traces[i].getHead().getFrame();
                break;
            case // Sort by size descending
            3:
                values[i] = -traces[i].size();
                break;
            case // Sort by length descending
            4:
                values[i] = -roi.getLength();
                break;
            case // Mean Square Displacement
            5:
                values[i] = -traces[i].getMSD();
                break;
            case // Mean / Frame
            6:
                values[i] = -traces[i].getMeanPerFrame();
                break;
        }
    }
    if (sort > 0)
        Sort.sort(indices, values);
    // Draw the traces as ROIs on an overlay
    Overlay o = new Overlay();
    LUT lut = LUTHelper.createLUT(DrawClusters.lut);
    final double scale = 256.0 / count;
    if (isUseStackPosition) {
        // Add the tracks on the frames containing the results
        final boolean isHyperStack = imp.isDisplayedHyperStack();
        for (int i = 0; i < count; i++) {
            final int index = indices[i];
            final Color c = LUTHelper.getColour(lut, (int) (i * scale));
            final PolygonRoi roi = (PolygonRoi) rois[index];
            roi.setFillColor(c);
            roi.setStrokeColor(c);
            final FloatPolygon fp = roi.getNonSplineFloatPolygon();
            // For each frame in the track, add the ROI track and a point ROI for the current position
            for (int j = 0; j < frames[index].length; j++) {
                addToOverlay(o, (Roi) roi.clone(), isHyperStack, frames[index][j]);
                //PointRoi pointRoi = new PointRoi(pos.x + fp.xpoints[j], pos.y + fp.ypoints[j]);
                PointRoi pointRoi = new PointRoi(fp.xpoints[j], fp.ypoints[j]);
                pointRoi.setPointType(3);
                pointRoi.setFillColor(c);
                pointRoi.setStrokeColor(Color.black);
                addToOverlay(o, pointRoi, isHyperStack, frames[index][j]);
            }
        }
    } else {
        // Add the tracks as a single overlay
        for (int i = 0; i < count; i++) {
            final Roi roi = rois[indices[i]];
            roi.setStrokeColor(new Color(lut.getRGB((int) (i * scale))));
            o.add(roi);
        }
    }
    imp.setOverlay(o);
    IJ.showStatus(msg);
}

Example 3 with ImageStack

use of ij.ImageStack in project GDSC-SMLM by aherbert.

the class SpotAnalysis method createProfile.

private void createProfile(ImagePlus imp, Rectangle bounds, double psfWidth, double blur) {
    areaBounds = bounds;
    this.imp = imp;
    area = bounds.width * bounds.height;
    clearSelectedFrames();
    // Get a profile through the images
    IJ.showStatus("Calculating raw profile");
    final int nSlices = imp.getStackSize();
    ImageStack rawSpot = new ImageStack(bounds.width, bounds.height, nSlices);
    double[][] profile = extractSpotProfile(imp, bounds, rawSpot);
    // Retain the existing display range
    double min = 0, max = Double.POSITIVE_INFINITY;
    if (rawImp != null) {
        min = rawImp.getDisplayRangeMin();
        max = rawImp.getDisplayRangeMax();
    }
    rawImp = showSpot(rawSpotTitle, rawSpot);
    if (max != Double.POSITIVE_INFINITY) {
        rawImp.setDisplayRange(min, max);
    }
    rawMean = profile[0];
    rawSd = profile[1];
    // Check if there are fitted results in memory
    addCandidateFrames(imp.getTitle());
    updateProfilePlots();
    if (blur > 0) {
        IJ.showStatus("Calculating blur ...");
        ImageStack stack = imp.getImageStack();
        ImageStack newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
        // Multi-thread the blur stage
        ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
        List<Future<?>> futures = new LinkedList<Future<?>>();
        Utils.setShowProgress(false);
        blurCount = 0;
        // TODO - See if this is faster if processing multiple slices in each worker
        int slices = 5;
        for (int n = 1; n <= nSlices; n += slices) {
            futures.add(threadPool.submit(new BlurWorker(stack, n, slices, bounds, blur * psfWidth, newStack)));
        }
        IJ.showStatus("Calculating blur ... Finishing");
        Utils.waitForCompletion(futures);
        threadPool.shutdown();
        Utils.setShowProgress(false);
        IJ.showStatus("Calculating blur ... Drawing");
        ImageStack blurSpot = new ImageStack(bounds.width, bounds.height, nSlices);
        extractSpotProfile(new ImagePlus("Blur", newStack), bounds, blurSpot);
        // Retain the existing display range
        max = Double.POSITIVE_INFINITY;
        if (blurImp != null) {
            min = blurImp.getDisplayRangeMin();
            max = blurImp.getDisplayRangeMax();
        }
        blurImp = showSpot(blurSpotTitle, blurSpot);
        if (max != Double.POSITIVE_INFINITY) {
            blurImp.setDisplayRange(min, max);
        }
        IJ.showStatus("");
    } else {
        blurImp = null;
    }
    // Add a z-projection of the blur/original image
    ZProjector project = new ZProjector((blurImp == null) ? rawImp : blurImp);
    project.setMethod(ZProjector.AVG_METHOD);
    project.doProjection();
    showSpot(avgSpotTitle, project.getProjection().getImageStack());
    if (!candidateFrames.isEmpty())
        // Set the first candidate frame
        rawImp.setSlice(candidateFrames.get(0));
    else
        updateCurrentSlice(rawImp.getCurrentSlice());
    IJ.showStatus("");
}

Example 4 with ImageStack

use of ij.ImageStack in project GDSC-SMLM by aherbert.

the class CreateData method drawImage.

//StoredDataStatistics rawPhotons = new StoredDataStatistics();
//StoredDataStatistics drawPhotons = new StoredDataStatistics();
//	private synchronized void addRaw(double d)
//	{
//		//rawPhotons.add(d);
//	}
//
//	private synchronized void addDraw(double d)
//	{
//		//drawPhotons.add(d);
//	}
/**
	 * Create an image from the localisations using the configured PSF width. Draws a new stack
	 * image.
	 * <p>
	 * Note that the localisations are filtered using the signal. The input list of localisations will be updated.
	 * 
	 * @param localisationSets
	 * @return The localisations
	 */
private List<LocalisationModel> drawImage(final List<LocalisationModelSet> localisationSets) {
    if (localisationSets.isEmpty())
        return null;
    // Create a new list for all localisation that are drawn (i.e. pass the signal filters)
    List<LocalisationModelSet> newLocalisations = Collections.synchronizedList(new ArrayList<LocalisationModelSet>(localisationSets.size()));
    photonsRemoved = new AtomicInteger();
    t1Removed = new AtomicInteger();
    tNRemoved = new AtomicInteger();
    photonStats = new SummaryStatistics();
    // Add drawn spots to memory
    results = new MemoryPeakResults();
    Calibration c = new Calibration(settings.pixelPitch, settings.getTotalGain(), settings.exposureTime);
    c.setEmCCD((settings.getEmGain() > 1));
    c.setBias(settings.bias);
    c.setReadNoise(settings.readNoise * ((settings.getCameraGain() > 0) ? settings.getCameraGain() : 1));
    c.setAmplification(settings.getAmplification());
    results.setCalibration(c);
    results.setSortAfterEnd(true);
    results.begin();
    maxT = localisationSets.get(localisationSets.size() - 1).getTime();
    // Display image
    ImageStack stack = new ImageStack(settings.size, settings.size, maxT);
    final double psfSD = getPsfSD();
    if (psfSD <= 0)
        return null;
    ImagePSFModel imagePSFModel = null;
    if (imagePSF) {
        // Create one Image PSF model that can be copied
        imagePSFModel = createImagePSF(localisationSets);
        if (imagePSFModel == null)
            return null;
    }
    IJ.showStatus("Drawing image ...");
    // Multi-thread for speed
    // Note that the default Executors.newCachedThreadPool() will continue to make threads if
    // new tasks are added. We need to limit the tasks that can be added using a fixed size
    // blocking queue.
    // http://stackoverflow.com/questions/1800317/impossible-to-make-a-cached-thread-pool-with-a-size-limit
    // ExecutorService threadPool = Executors.newCachedThreadPool();
    ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    List<Future<?>> futures = new LinkedList<Future<?>>();
    // Count all the frames to process
    frame = 0;
    totalFrames = maxT;
    // Collect statistics on the number of photons actually simulated
    // Process all frames
    int i = 0;
    int lastT = -1;
    for (LocalisationModelSet l : localisationSets) {
        if (Utils.isInterrupted())
            break;
        if (l.getTime() != lastT) {
            lastT = l.getTime();
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, i, lastT, createPSFModel(imagePSFModel), results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
        i++;
    }
    // Finish processing data
    Utils.waitForCompletion(futures);
    futures.clear();
    if (Utils.isInterrupted()) {
        IJ.showProgress(1);
        return null;
    }
    // Do all the frames that had no localisations
    for (int t = 1; t <= maxT; t++) {
        if (Utils.isInterrupted())
            break;
        if (stack.getPixels(t) == null) {
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, maxT, t, null, results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
    }
    // Finish
    Utils.waitForCompletion(futures);
    threadPool.shutdown();
    IJ.showProgress(1);
    if (Utils.isInterrupted()) {
        return null;
    }
    results.end();
    // Clear memory
    imagePSFModel = null;
    threadPool = null;
    futures.clear();
    futures = null;
    if (photonsRemoved.get() > 0)
        Utils.log("Removed %d localisations with less than %.1f rendered photons", photonsRemoved.get(), settings.minPhotons);
    if (t1Removed.get() > 0)
        Utils.log("Removed %d localisations with no neighbours @ SNR %.2f", t1Removed.get(), settings.minSNRt1);
    if (tNRemoved.get() > 0)
        Utils.log("Removed %d localisations with valid neighbours @ SNR %.2f", tNRemoved.get(), settings.minSNRtN);
    if (photonStats.getN() > 0)
        Utils.log("Average photons rendered = %s +/- %s", Utils.rounded(photonStats.getMean()), Utils.rounded(photonStats.getStandardDeviation()));
    //System.out.printf("rawPhotons = %f\n", rawPhotons.getMean());
    //System.out.printf("drawPhotons = %f\n", drawPhotons.getMean());
    //Utils.showHistogram("draw photons", drawPhotons, "photons", true, 0, 1000);
    // Update with all those localisation that have been drawn
    localisationSets.clear();
    localisationSets.addAll(newLocalisations);
    newLocalisations = null;
    IJ.showStatus("Displaying image ...");
    ImageStack newStack = stack;
    if (!settings.rawImage) {
        // Get the global limits and ensure all values can be represented
        Object[] imageArray = stack.getImageArray();
        float[] limits = Maths.limits((float[]) imageArray[0]);
        for (int j = 1; j < imageArray.length; j++) limits = Maths.limits(limits, (float[]) imageArray[j]);
        // Leave bias in place
        limits[0] = 0;
        // Check if the image will fit in a 16-bit range
        if ((limits[1] - limits[0]) < 65535) {
            // Convert to 16-bit
            newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
            // Account for rounding
            final float min = (float) (limits[0] - 0.5);
            for (int j = 0; j < imageArray.length; j++) {
                float[] image = (float[]) imageArray[j];
                short[] pixels = new short[image.length];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = (short) (image[k] - min);
                }
                newStack.setPixels(pixels, j + 1);
                // Free memory
                imageArray[j] = null;
                // Attempt to stay within memory (check vs 32MB)
                if (MemoryPeakResults.freeMemory() < 33554432L)
                    MemoryPeakResults.runGCOnce();
            }
        } else {
            // Keep as 32-bit but round to whole numbers
            for (int j = 0; j < imageArray.length; j++) {
                float[] pixels = (float[]) imageArray[j];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = Math.round(pixels[k]);
                }
            }
        }
    }
    // Show image
    ImagePlus imp = Utils.display(CREATE_DATA_IMAGE_TITLE, newStack);
    ij.measure.Calibration cal = new ij.measure.Calibration();
    String unit = "nm";
    double unitPerPixel = settings.pixelPitch;
    if (unitPerPixel > 100) {
        unit = "um";
        unitPerPixel /= 1000.0;
    }
    cal.setUnit(unit);
    cal.pixelHeight = cal.pixelWidth = unitPerPixel;
    imp.setCalibration(cal);
    imp.setDimensions(1, 1, newStack.getSize());
    imp.resetDisplayRange();
    imp.updateAndDraw();
    saveImage(imp);
    results.setSource(new IJImageSource(imp));
    results.setName(CREATE_DATA_IMAGE_TITLE + " (" + TITLE + ")");
    results.setConfiguration(createConfiguration((float) psfSD));
    results.setBounds(new Rectangle(0, 0, settings.size, settings.size));
    MemoryPeakResults.addResults(results);
    setBenchmarkResults(imp, results);
    if (benchmarkMode && benchmarkParameters != null)
        benchmarkParameters.setPhotons(results);
    List<LocalisationModel> localisations = toLocalisations(localisationSets);
    savePulses(localisations, results, CREATE_DATA_IMAGE_TITLE);
    // Saved the fixed and moving localisations into different datasets
    saveFixedAndMoving(results, CREATE_DATA_IMAGE_TITLE);
    return localisations;
}

Example 5 with ImageStack

use of ij.ImageStack in project GDSC-SMLM by aherbert.

the class CreateData method createMaskDistribution.

private SpatialDistribution createMaskDistribution(ImagePlus imp, double sliceDepth, boolean updateArea) {
    // Calculate the scale of the mask
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    final double scaleX = (double) settings.size / w;
    final double scaleY = (double) settings.size / h;
    // Use an image for the distribution
    if (imp.getStackSize() > 1) {
        ImageStack stack = imp.getImageStack();
        List<int[]> masks = new ArrayList<int[]>(stack.getSize());
        int[] maxMask = new int[w * h];
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            int[] mask = extractMask(stack.getProcessor(slice));
            if (updateArea) {
                for (int i = 0; i < mask.length; i++) if (mask[i] != 0)
                    maxMask[i] = 1;
            }
            masks.add(mask);
        }
        if (updateArea)
            updateArea(maxMask, w, h);
        if (sliceDepth == 0)
            // Auto configure to the full depth of the simulation
            sliceDepth = settings.depth / masks.size();
        return new MaskDistribution3D(masks, w, h, sliceDepth / settings.pixelPitch, scaleX, scaleY, createRandomGenerator());
    } else {
        int[] mask = extractMask(imp.getProcessor());
        if (updateArea)
            updateArea(mask, w, h);
        return new MaskDistribution(mask, w, h, settings.depth / settings.pixelPitch, scaleX, scaleY, createRandomGenerator());
    }
}