Examples with FitMode uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FitMode used on opensource projects

Example 1 with FitMode

use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.FitMode 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;
}