Search in sources :

Example 1 with Spot

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot in project GDSC-SMLM by aherbert.

the class TsfPeakResultsWriter method add.

@Override
public void add(int peak, int origX, int origY, float origValue, double error, float noise, float meanIntensity, float[] params, float[] paramsStdDev) {
    if (out == null) {
        return;
    }
    final Spot.Builder builder = Spot.newBuilder();
    builder.setMolecule(id.incrementAndGet());
    builder.setChannel(1);
    builder.setFluorophoreType(1);
    builder.setFrame(peak);
    builder.setXPosition(origX);
    builder.setYPosition(origY);
    setParam(params, builder);
    builder.setError(error);
    builder.setNoise(noise);
    builder.setMeanIntensity(meanIntensity);
    builder.setOriginalValue(origValue);
    if (paramsStdDev != null) {
        addNewParamStdDevs(builder, paramsStdDev);
    }
    final Spot spot = builder.build();
    writeResult(1, spot);
}
Also used : Spot(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot) Builder(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot.Builder)

Example 2 with Spot

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot in project GDSC-SMLM by aherbert.

the class ResultsManagerTest method checkWriteTsfMatchesRead.

@SuppressWarnings("null")
private static void checkWriteTsfMatchesRead(RandomSeed seed, int channels, int slices, int positions, int types) {
    Assumptions.assumeFalse(java.awt.GraphicsEnvironment.isHeadless());
    final String filename = createFile();
    FileOutputStream out = null;
    try {
        out = new FileOutputStream(filename);
    } catch (final Exception ex) {
        closeOutput(out);
        ex.printStackTrace();
        Assertions.fail(ex.getMessage());
    }
    // Write the offsets used in the Tsf format
    try {
        @SuppressWarnings("resource") final DataOutputStream dos = new DataOutputStream(out);
        dos.writeInt(0);
        dos.writeLong(0);
    } catch (final IOException ex) {
        closeOutput(out);
        ex.printStackTrace();
        Assertions.fail(ex.getMessage());
    }
    // Generate random spots
    final UniformRandomProvider rand = RngUtils.create(seed.getSeed());
    final int size = 100;
    final Spot[] spots = new Spot[size];
    for (int i = 1; i <= size; i++) {
        final Spot.Builder builder = Spot.newBuilder();
        builder.setChannel(nextInt(rand, channels));
        builder.setSlice(nextInt(rand, slices));
        builder.setPos(nextInt(rand, positions));
        builder.setFluorophoreType(nextInt(rand, types));
        // This is a required field but is ignored when reading
        builder.setMolecule(i);
        builder.setCluster(rand.nextInt(10));
        builder.setFrame(nextInt(rand, 100));
        builder.setXPosition(rand.nextInt(50));
        builder.setYPosition(rand.nextInt(50));
        builder.setBackground(rand.nextFloat());
        builder.setIntensity(rand.nextFloat());
        builder.setX(rand.nextFloat());
        builder.setY(rand.nextFloat());
        builder.setZ(rand.nextFloat());
        builder.setWidth((float) (Gaussian2DFunction.SD_TO_FWHM_FACTOR * rand.nextDouble()));
        final Spot spot = builder.build();
        spots[i - 1] = spot;
        try {
            spot.writeDelimitedTo(out);
        } catch (final IOException ex) {
            closeOutput(out);
            ex.printStackTrace();
            Assertions.fail(ex.getMessage());
        }
    }
    // Write the header
    // Get the offset to the SpotList message
    long offset = 0;
    try {
        // The offset is the amount to skip forward after reading the int
        // magic number (4 bytes) and long offset (8 bytes)
        // out.flush();
        offset = out.getChannel().position() - 12;
    } catch (final IOException ex) {
        closeOutput(out);
        ex.printStackTrace();
        Assertions.fail(ex.getMessage());
    }
    // Record the SpotList message
    final SpotList.Builder builder = SpotList.newBuilder();
    builder.setApplicationId(1);
    builder.setNrSpots(size);
    builder.setLocationUnits(LocationUnits.PIXELS);
    builder.setIntensityUnits(IntensityUnits.COUNTS);
    builder.setFitMode(FitMode.ONEAXIS);
    builder.setNrChannels(channels);
    builder.setNrSlices(slices);
    builder.setNrPos(positions);
    for (int type = 1; type <= types; type++) {
        final FluorophoreType.Builder typeBuilder = FluorophoreType.newBuilder();
        typeBuilder.setId(type);
        typeBuilder.setDescription("Type " + type);
        typeBuilder.setIsFiducial(rand.nextDouble() < 0.5);
        builder.addFluorophoreTypes(typeBuilder.build());
    }
    final SpotList spotList = builder.build();
    try {
        spotList.writeDelimitedTo(out);
    } catch (final IOException ex) {
        ex.printStackTrace();
        Assertions.fail(ex.getMessage());
    } finally {
        closeOutput(out);
    }
    // Write the offset to the SpotList message into the offset position
    try (RandomAccessFile f = new RandomAccessFile(new File(filename), "rw")) {
        f.seek(4);
        f.writeLong(offset);
    } catch (final Exception ex) {
        ex.printStackTrace();
        Assertions.fail(ex.getMessage());
    }
    // Read each combination
    for (int channel = 1; channel <= channels; channel++) {
        for (int slice = 1; slice <= slices; slice++) {
            for (int position = 1; position <= positions; position++) {
                for (int type = 1; type <= types; type++) {
                    final StringBuilder sb = new StringBuilder();
                    sb.append(" channel=").append(channel);
                    sb.append(" slice=").append(slice);
                    sb.append(" position=").append(position);
                    sb.append(" fluorophore_type=[").append(type).append(":Type ").append(type).append(":fiducial=").append(builder.getFluorophoreTypes(type - 1).getIsFiducial()).append(']');
                    // This is needed to trick the Macro class into returning the options
                    // for the thread to the GenericDialog used in the ResultsManager
                    Thread.currentThread().setName("Run$_");
                    Macro.setOptions(sb.toString());
                    final MemoryPeakResults in = ResultsManager.loadInputResults(ResultsManager.INPUT_FILE, false, null, null, new ResultsManager.FilenameLoadOption(filename));
                    checkEqual(spots, channel, slice, position, type, in);
                }
            }
        }
    }
}
Also used : FluorophoreType(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FluorophoreType) Spot(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot) DataOutputStream(java.io.DataOutputStream) IOException(java.io.IOException) SpotList(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.SpotList) IOException(java.io.IOException) RandomAccessFile(java.io.RandomAccessFile) FileOutputStream(java.io.FileOutputStream) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults) RandomAccessFile(java.io.RandomAccessFile) File(java.io.File)

Example 3 with Spot

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot in project GDSC-SMLM by aherbert.

the class ResultsManagerTest method extract.

private static MemoryPeakResults extract(Spot[] spots, int channel, int slice, int position, int type) {
    final MemoryPeakResults results = new MemoryPeakResults(PsfHelper.create(PSFType.ONE_AXIS_GAUSSIAN_2D));
    for (final Spot spot : spots) {
        if (spot.getChannel() == channel && spot.getSlice() == slice && spot.getPos() == position && spot.getFluorophoreType() == type) {
            final int id = spot.getCluster();
            final int startFrame = spot.getFrame();
            final int origX = spot.getXPosition();
            final int origY = spot.getYPosition();
            final float origValue = 0;
            final double error = 0;
            final float noise = 0;
            final float[] params = Gaussian2DPeakResultHelper.createOneAxisParams(spot.getBackground(), spot.getIntensity(), spot.getX(), spot.getY(), spot.getZ(), (float) (spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR));
            final float[] paramsStdDev = null;
            final IdPeakResult peak = new IdPeakResult(startFrame, origX, origY, origValue, error, noise, 0, params, paramsStdDev, id);
            results.add(peak);
        }
    }
    return results;
}
Also used : Spot(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot) IdPeakResult(uk.ac.sussex.gdsc.smlm.results.IdPeakResult) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)

Example 4 with Spot

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot in project GDSC-SMLM by aherbert.

the class TsfPeakResultsReader method read.

/**
 * Read the results from the TSF file into memory.
 *
 * @return The results set (or null if an error occurred)
 */
public MemoryPeakResults read() {
    readHeader();
    if (spotList == null) {
        return null;
    }
    final MemoryPeakResults results = createResults();
    // Used in the exception handler to check the correct number of spots were read
    long expectedSpots = -1;
    // Read the messages that contain the spot data
    try (BufferedInputStream fi = new BufferedInputStream(new FileInputStream(filename))) {
        // size of int + size of long
        FileUtils.skip(fi, 12);
        final LocationUnits locationUnits = spotList.getLocationUnits();
        boolean locationUnitsWarning = false;
        final IntensityUnits intensityUnits = spotList.getIntensityUnits();
        boolean intensityUnitsWarning = false;
        final FitMode fitMode = (spotList.hasFitMode()) ? spotList.getFitMode() : FitMode.ONEAXIS;
        final boolean filterPosition = position > 0;
        final boolean filterSlice = slice > 0;
        // Set up to read two-axis and theta data
        final PSF psf = PsfHelper.create(PSFType.TWO_AXIS_AND_THETA_GAUSSIAN_2D);
        final int[] indices = PsfHelper.getGaussian2DWxWyIndices(psf);
        final int isx = indices[0];
        final int isy = indices[1];
        final int ia = PsfHelper.getGaussian2DAngleIndex(psf);
        int nparams = PeakResult.STANDARD_PARAMETERS;
        switch(fitMode) {
            case ONEAXIS:
                nparams += 1;
                break;
            case TWOAXIS:
                nparams += 2;
                break;
            case TWOAXISANDTHETA:
                nparams += 3;
                break;
            default:
                logger.log(Level.WARNING, () -> "Unknown fit mode: " + fitMode);
                return null;
        }
        expectedSpots = getExpectedSpots();
        long totalSpots = 0;
        while (fi.available() > 0 && (totalSpots < expectedSpots || expectedSpots == 0)) {
            totalSpots++;
            final Spot spot = Spot.parseDelimitedFrom(fi);
            // Only read the specified channel, position, slice and fluorophore type
            if (spot.getChannel() != channel) {
                continue;
            }
            if (filterPosition && spot.hasPos() && spot.getPos() != position) {
                continue;
            }
            if (filterSlice && spot.hasSlice() && spot.getSlice() != slice) {
                continue;
            }
            if (spot.hasFluorophoreType() && spot.getFluorophoreType() != fluorophoreType) {
                continue;
            }
            if (spot.hasLocationUnits() && !locationUnitsWarning && spot.getLocationUnits() != locationUnits) {
                logger.warning(() -> "Spot message has different location units, the units will be ignored: " + spot.getLocationUnits());
                locationUnitsWarning = true;
            }
            if (spot.hasIntensityUnits() && !intensityUnitsWarning && spot.getIntensityUnits() != intensityUnits) {
                logger.warning(() -> "Spot message has different intensity units, the units will be ignored: " + spot.getIntensityUnits());
                intensityUnitsWarning = true;
            }
            // Required fields
            final int frame = spot.getFrame();
            final float[] params = new float[nparams];
            params[PeakResult.X] = spot.getX();
            params[PeakResult.Y] = spot.getY();
            params[PeakResult.INTENSITY] = spot.getIntensity();
            if (spot.hasBackground()) {
                params[PeakResult.BACKGROUND] = spot.getBackground();
            }
            if (spot.hasZ()) {
                params[PeakResult.Z] = spot.getZ();
            }
            // Support different Gaussian shapes
            if (fitMode == FitMode.ONEAXIS) {
                params[isx] = (float) (spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR);
            } else {
                if (!spot.hasA()) {
                    params[isx] = params[isy] = (float) (spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR);
                } else {
                    final double a = Math.sqrt(spot.getA());
                    final double sd = spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR;
                    params[isx] = (float) (sd * a);
                    params[isy] = (float) (sd / a);
                }
                if (fitMode == FitMode.TWOAXISANDTHETA && spot.hasTheta()) {
                    params[ia] = spot.getTheta();
                }
            }
            // We can use the original position in pixels used for fitting
            final int origX = (spot.hasXPosition()) ? spot.getXPosition() : 0;
            final int origY = (spot.hasYPosition()) ? spot.getYPosition() : 0;
            // Q. Should we use the required field 'molecule'?
            float origValue = 0;
            double error = 0;
            float noise = 0;
            float meanIntensity = 0;
            float[] paramsStdDev = null;
            int id = Integer.MIN_VALUE;
            int category = Integer.MIN_VALUE;
            int endFrame = Integer.MIN_VALUE;
            if (isGdsc) {
                error = spot.getError();
                noise = spot.getNoise();
                meanIntensity = spot.getMeanIntensity();
                id = spot.getCluster();
                category = spot.getCategory();
                origValue = spot.getOriginalValue();
                endFrame = spot.getEndFrame();
                if (spot.getParamStdDevsCount() != 0) {
                    paramsStdDev = new float[spot.getParamStdDevsCount()];
                    for (int i = 0; i < paramsStdDev.length; i++) {
                        paramsStdDev[i] = spot.getParamStdDevs(i);
                    }
                }
            } else if (spot.hasCluster()) {
                // Use the standard cluster field for the ID.
                // Note: The GDSC SMLM results do not distinguish molecule or cluster so here
                // we pick cluster first.
                id = spot.getCluster();
            } else if (spot.hasMolecule()) {
                // If no cluster then use the molecule.
                // Note: This may just be a natural series with a unique number for each localisation.
                id = spot.getMolecule();
            }
            // Allow storing any of the optional attributes
            final AttributePeakResult peakResult = new AttributePeakResult(frame, origX, origY, origValue, error, noise, meanIntensity, params, paramsStdDev);
            if (id != Integer.MIN_VALUE) {
                peakResult.setId(id);
            }
            if (category != Integer.MIN_VALUE) {
                peakResult.setCategory(category);
            }
            if (endFrame != Integer.MIN_VALUE) {
                peakResult.setEndFrame(endFrame);
            }
            if (spot.hasXPrecision() || spot.hasYPrecision()) {
                // Use the average. Note this is not the Euclidean distance since we divide by n
                double sumSq = 0;
                int count = 0;
                if (spot.hasXPrecision()) {
                    sumSq = spot.getXPrecision() * spot.getXPrecision();
                    count++;
                }
                if (spot.hasYPrecision()) {
                    sumSq += spot.getYPrecision() * spot.getYPrecision();
                    count++;
                }
                peakResult.setPrecision(Math.sqrt(sumSq / count));
            }
            results.add(peakResult);
        }
    } catch (final IOException ex) {
        logger.log(Level.WARNING, ex, () -> "Failed to read TSF file: " + filename);
        if (expectedSpots == -1) {
            // The exception was created during set-up.
            return null;
        }
        // Only fail if there is a number of expected spots.
        if (expectedSpots != 0) {
            logger.warning(() -> "Unexpected error in reading Spot messages, no results will be returned");
            return null;
        }
    }
    return results;
}
Also used : PSF(uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF) LocationUnits(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.LocationUnits) Spot(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot) IntensityUnits(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.IntensityUnits) IOException(java.io.IOException) FileInputStream(java.io.FileInputStream) BufferedInputStream(java.io.BufferedInputStream) FitMode(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FitMode)

Example 5 with Spot

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot in project GDSC-SMLM by aherbert.

the class TsfPeakResultsWriter method addAll.

@Override
public void addAll(PeakResult[] results) {
    if (out == null) {
        return;
    }
    final Spot[] spots = new Spot[20];
    int count = 0;
    final Spot.Builder builder = Spot.newBuilder();
    for (final PeakResult result : results) {
        addStandardFields(builder, result);
        addParamStdDevs(builder, result.getParameterDeviations());
        spots[count++] = builder.build();
        // Flush the output to allow for very large input lists
        if (count >= spots.length) {
            writeResult(count, spots);
            if (!isActive()) {
                return;
            }
            count = 0;
        }
    }
    writeResult(count, spots);
}
Also used : Spot(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot) Builder(uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot.Builder)

Aggregations

Spot (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot)5 IOException (java.io.IOException)2 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)2 Builder (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.Spot.Builder)2 BufferedInputStream (java.io.BufferedInputStream)1 DataOutputStream (java.io.DataOutputStream)1 File (java.io.File)1 FileInputStream (java.io.FileInputStream)1 FileOutputStream (java.io.FileOutputStream)1 RandomAccessFile (java.io.RandomAccessFile)1 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)1 PSF (uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF)1 IdPeakResult (uk.ac.sussex.gdsc.smlm.results.IdPeakResult)1 FitMode (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FitMode)1 FluorophoreType (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FluorophoreType)1 IntensityUnits (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.IntensityUnits)1 LocationUnits (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.LocationUnits)1 SpotList (uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.SpotList)1