Examples with PSF uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF used on opensource projects

Example 11 with PSF

use of uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF in project GDSC-SMLM by aherbert.

the class TsfPeakResultsReader method createResults.

private MemoryPeakResults createResults() {
    // Limit the capacity since we may not need all the spots
    int capacity = 1000;
    if (spotList.hasNrSpots()) {
        capacity = (int) Math.min(100000, spotList.getNrSpots());
    }
    final MemoryPeakResults results = new MemoryPeakResults(capacity);
    // Create the type of Gaussian PSF
    if (spotList.hasFitMode()) {
        switch(spotList.getFitMode()) {
            case ONEAXIS:
                results.setPsf(PsfHelper.create(PSFType.ONE_AXIS_GAUSSIAN_2D));
                break;
            case TWOAXIS:
                results.setPsf(PsfHelper.create(PSFType.TWO_AXIS_GAUSSIAN_2D));
                break;
            case TWOAXISANDTHETA:
                results.setPsf(PsfHelper.create(PSFType.TWO_AXIS_AND_THETA_GAUSSIAN_2D));
                break;
            default:
                break;
        }
    }
    // Generic reconstruction
    String name;
    if (spotList.hasName()) {
        name = spotList.getName();
    } else {
        name = FileUtils.getName(filename);
    }
    // Append these if not using the defaults
    if (channel != 1 || slice != 0 || position != 0 || fluorophoreType != 1) {
        name = String.format("%s c=%d, s=%d, p=%d, ft=%d", name, channel, slice, position, fluorophoreType);
    }
    results.setName(name);
    // if (spotList.hasNrPixelsX() && spotList.hasNrPixelsY())
    // {
    // // Do not do this. The size of the camera may not map to the data bounds due
    // // to the support for position offsets.
    // results.setBounds(new Rectangle(0, 0, spotList.getNrPixelsX(), spotList.getNrPixelsY()));
    // }
    final CalibrationWriter cal = new CalibrationWriter();
    // Spots are associated with frames
    cal.setTimeUnit(TimeUnit.FRAME);
    if (spotList.hasPixelSize()) {
        cal.setNmPerPixel(spotList.getPixelSize());
    }
    if (spotList.getEcfCount() >= channel) {
        // ECF is per channel
        final double ecf = spotList.getEcf(channel - 1);
        // QE is per fluorophore type
        final double qe = (spotList.getQeCount() >= fluorophoreType) ? spotList.getQe(fluorophoreType - 1) : 1;
        // e-/photon / e-/count => count/photon
        cal.setCountPerPhoton(qe / ecf);
        cal.setQuantumEfficiency(qe);
    }
    if (isGdsc) {
        if (spotList.hasSource()) {
            // Deserialise
            results.setSource(ImageSource.fromXml(spotList.getSource()));
        }
        if (spotList.hasRoi()) {
            final ROI roi = spotList.getRoi();
            if (roi.hasX() && roi.hasY() && roi.hasXWidth() && roi.hasYWidth()) {
                results.setBounds(new Rectangle(roi.getX(), roi.getY(), roi.getXWidth(), roi.getYWidth()));
            }
        }
        if (spotList.hasGain()) {
            cal.setCountPerPhoton(spotList.getGain());
        }
        if (spotList.hasExposureTime()) {
            cal.setExposureTime(spotList.getExposureTime());
        }
        if (spotList.hasReadNoise()) {
            cal.setReadNoise(spotList.getReadNoise());
        }
        if (spotList.hasBias()) {
            cal.setBias(spotList.getBias());
        }
        if (spotList.hasCameraType()) {
            cal.setCameraType(cameraTypeMap.get(spotList.getCameraType()));
        } else {
            cal.setCameraType(null);
        }
        if (spotList.hasConfiguration()) {
            results.setConfiguration(spotList.getConfiguration());
        }
        // Allow restoring the GDSC PSF exactly
        if (spotList.hasPSF()) {
            try {
                final Parser parser = JsonFormat.parser();
                final PSF.Builder psfBuilder = PSF.newBuilder();
                parser.merge(spotList.getPSF(), psfBuilder);
                results.setPsf(psfBuilder.build());
            } catch (final InvalidProtocolBufferException ex) {
                logger.warning("Unable to deserialise the PSF settings");
            }
        }
    }
    if (spotList.hasLocationUnits()) {
        cal.setDistanceUnit(locationUnitsMap.get(spotList.getLocationUnits()));
        if (!spotList.hasPixelSize() && spotList.getLocationUnits() != LocationUnits.PIXELS) {
            logger.warning(() -> "TSF location units are not pixels and no pixel size calibration is available." + " The dataset will be constructed in the native units: " + spotList.getLocationUnits());
        }
    } else {
        cal.setDistanceUnit(null);
    }
    if (spotList.hasIntensityUnits()) {
        cal.setIntensityUnit(intensityUnitsMap.get(spotList.getIntensityUnits()));
        if (!spotList.hasGain() && spotList.getIntensityUnits() != IntensityUnits.COUNTS) {
            logger.warning(() -> "TSF intensity units are not counts and no gain calibration is available." + " The dataset will be constructed in the native units: " + spotList.getIntensityUnits());
        }
    } else {
        cal.setIntensityUnit(null);
    }
    if (spotList.hasThetaUnits()) {
        cal.setAngleUnit(thetaUnitsMap.get(spotList.getThetaUnits()));
    } else {
        cal.setAngleUnit(null);
    }
    results.setCalibration(cal.getCalibration());
    return results;
}

Example 12 with PSF

use of uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF in project GDSC-SMLM by aherbert.

the class PsfProtosHelperTest method canConvertAstigmatismModel.

@Test
void canConvertAstigmatismModel() {
    // Use a reasonable z-depth function from the Smith, et al (2010) paper (page 377)
    final double sx = 1.08;
    final double sy = 1.01;
    final double gamma = 0.389;
    final double d = 0.531;
    final double Ax = -0.0708;
    final double Bx = -0.073;
    final double Ay = 0.164;
    final double By = 0.0417;
    final double nmPerPixel = 100;
    // Ax = Ay = 0;
    // Bx = By = 0;
    final DistanceUnit zDistanceUnit = DistanceUnit.UM;
    final DistanceUnit sDistanceUnit = DistanceUnit.PIXEL;
    final AstigmatismModel.Builder builder = AstigmatismModel.newBuilder();
    builder.setGamma(gamma);
    builder.setD(d);
    builder.setS0X(sx);
    builder.setAx(Ax);
    builder.setBx(Bx);
    builder.setS0Y(sy);
    builder.setAy(Ay);
    builder.setBy(By);
    builder.setZDistanceUnit(zDistanceUnit);
    builder.setSDistanceUnit(sDistanceUnit);
    builder.setNmPerPixel(nmPerPixel);
    final AstigmatismModel model1 = builder.build();
    final PSF psf = PsfProtosHelper.createPsf(model1, zDistanceUnit, sDistanceUnit);
    final AstigmatismModel model2 = PsfProtosHelper.createModel(psf, zDistanceUnit, sDistanceUnit, nmPerPixel);
    Assertions.assertEquals(model1, model2);
}

Example 13 with PSF

use of uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF in project GDSC-SMLM by aherbert.

the class SpotFinderPreview method run.

@Override
public void run(ImageProcessor ip) {
    if (refreshing) {
        return;
    }
    final Rectangle bounds = ip.getRoi();
    // Only do this if the settings changed
    final Calibration calibration = fitConfig.getCalibration();
    final FitEngineSettings fitEngineSettings = config.getFitEngineSettings();
    final PSF psf = fitConfig.getPsf();
    boolean newCameraModel = filter == null;
    if (!calibration.equals(lastCalibration)) {
        newCameraModel = true;
        // Set a camera model.
        // We have to set the camera type too to avoid configuration errors.
        CameraModel cameraModel = CameraModelManager.load(fitConfig.getCameraModelName());
        if (cameraModel == null) {
            cameraModel = new FakePerPixelCameraModel(0, 1, 1);
            fitConfig.setCameraType(CameraType.EMCCD);
        } else {
            fitConfig.setCameraType(CameraType.SCMOS);
            // Support cropped origin selection.
            final Rectangle sourceBounds = IJImageSource.getBounds(imp);
            cameraModel = PeakFit.cropCameraModel(cameraModel, sourceBounds, null, true);
            if (cameraModel == null) {
                gd.getPreviewCheckbox().setState(false);
                return;
            }
        }
        fitConfig.setCameraModel(cameraModel);
    }
    if (newCameraModel || !fitEngineSettings.equals(lastFitEngineSettings) || !psf.equals(lastPsf)) {
        // Configure a jury filter
        if (config.getDataFilterType() == DataFilterType.JURY && !PeakFit.configureDataFilter(config, PeakFit.FLAG_NO_SAVE)) {
            gd.getPreviewCheckbox().setState(false);
            return;
        }
        try {
            filter = config.createSpotFilter();
        } catch (final Exception ex) {
            filter = null;
            this.imp.setOverlay(overlay);
            // Required for ImageJ to disable the preview
            throw new IllegalStateException("Unable to create spot filter", ex);
        }
        ImageJUtils.log(filter.getDescription());
    }
    lastCalibration = calibration;
    lastFitEngineSettings = fitEngineSettings;
    lastPsf = psf;
    // This code can probably be removed since the crop is done above.
    if (fitConfig.getCameraTypeValue() == CameraType.SCMOS_VALUE) {
        // Instead just warn if the roi cannot be extracted from the selected model
        // or there is a mismatch
        final Rectangle modelBounds = fitConfig.getCameraModel().getBounds();
        if (modelBounds != null) {
            if (!modelBounds.contains(bounds)) {
                // @formatter:off
                ImageJUtils.log("WARNING: Camera model bounds [x=%d,y=%d,width=%d,height=%d]" + " does not contain image target bounds [x=%d,y=%d,width=%d,height=%d]", modelBounds.x, modelBounds.y, modelBounds.width, modelBounds.height, bounds.x, bounds.y, bounds.width, bounds.height);
            // @formatter:on
            // Warn if the model bounds are mismatched than the image as this may be an incorrect
            // selection for the camera model
            } else if (modelBounds.x != 0 || modelBounds.y != 0 || modelBounds.width > ip.getWidth() || modelBounds.height > ip.getHeight()) {
                // @formatter:off
                ImageJUtils.log("WARNING: Probably an incorrect camera model!\n" + "Model bounds [x=%d,y=%d,width=%d,height=%d]\n" + "do not match the image target bounds [width=%d,height=%d].", modelBounds.x, modelBounds.y, modelBounds.width, modelBounds.height, ip.getWidth(), ip.getHeight());
            // @formatter:on
            }
        }
    }
    run(ip, filter);
}

Example 14 with PSF

use of uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF in project GDSC-SMLM by aherbert.

the class SpotInspector method getStandardDeviation.

private float getStandardDeviation(MemoryPeakResults results2) {
    // Standard deviation is only needed for the width filtering
    if (settings.sortOrderIndex != 8) {
        return 0;
    }
    final PSF psf = results2.getPsf();
    if (!PsfHelper.isGaussian2D(psf)) {
        return -1;
    }
    final FitConfiguration fitConfig = new FitConfiguration();
    fitConfig.setPsf(psf);
    return (float) MathUtils.max(1, fitConfig.getInitialXSd(), fitConfig.getInitialYSd());
}

Example 15 with PSF

use of uk.ac.sussex.gdsc.smlm.data.config.PSFProtos.PSF 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;
}