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

Example 1 with CalibrationReader

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

the class DoubletAnalysis method summariseResults.

/**
 * Summarise results.
 *
 * @param results the results
 * @param density the density
 * @param runTime the run time
 */
private void summariseResults(ArrayList<DoubletResult> results, double density, long runTime) {
    // If we are only assessing results with no neighbour candidates
    // TODO - Count the number of actual results that have no neighbours
    final int numberOfMolecules = this.results.size() - ignored.get();
    final FitConfiguration fitConfig = config.getFitConfiguration();
    // Store details we want in the analysis table
    final StringBuilder sb = new StringBuilder();
    sb.append(MathUtils.rounded(density)).append('\t');
    sb.append(MathUtils.rounded(getSa())).append('\t');
    sb.append(config.getFittingWidth()).append('\t');
    sb.append(PsfProtosHelper.getName(fitConfig.getPsfType()));
    sb.append(":").append(PeakFit.getSolverName(fitConfig));
    if (fitConfig.isModelCameraMle()) {
        sb.append(":Camera\t");
        // Add details of the noise model for the MLE
        final CalibrationReader r = new CalibrationReader(fitConfig.getCalibration());
        sb.append("EM=").append(r.isEmCcd());
        sb.append(":A=").append(MathUtils.rounded(r.getCountPerElectron()));
        sb.append(":N=").append(MathUtils.rounded(r.getReadNoise()));
        sb.append('\t');
    } else {
        sb.append("\t\t");
    }
    final String analysisPrefix = sb.toString();
    // -=-=-=-=-
    showResults(results, settings.showResults);
    sb.setLength(0);
    final int n = countN(results);
    // Create the benchmark settings and the fitting settings
    sb.append(numberOfMolecules).append('\t');
    sb.append(n).append('\t');
    sb.append(MathUtils.rounded(density)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.minSignal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.maxSignal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.averageSignal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.pixelPitch)).append('\t');
    sb.append(MathUtils.rounded(getSa() * simulationParameters.pixelPitch)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.gain)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.readNoise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.background)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.noise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.averageSignal / simulationParameters.noise)).append('\t');
    sb.append(config.getFittingWidth()).append('\t');
    sb.append(PsfProtosHelper.getName(fitConfig.getPsfType()));
    sb.append(":").append(PeakFit.getSolverName(fitConfig));
    if (fitConfig.isModelCameraMle()) {
        sb.append(":Camera\t");
        // Add details of the noise model for the MLE
        final CalibrationReader r = new CalibrationReader(fitConfig.getCalibration());
        sb.append("EM=").append(r.isEmCcd());
        sb.append(":A=").append(MathUtils.rounded(r.getCountPerElectron()));
        sb.append(":N=").append(MathUtils.rounded(r.getReadNoise()));
        sb.append('\t');
    } else {
        sb.append("\t\t");
    }
    // Now output the actual results ...
    // Show histograms as cumulative to avoid problems with bin width
    // Residuals scores
    // Iterations and evaluations where fit was OK
    final StoredDataStatistics[] stats = new StoredDataStatistics[Settings.NAMES2.length];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = new StoredDataStatistics();
    }
    // For Jaccard scoring we need to count the score with no residuals threshold,
    // i.e. Accumulate the score accepting all doublets that were fit
    double tp = 0;
    double fp = 0;
    double bestTp = 0;
    double bestFp = 0;
    final ArrayList<DoubletBonus> data = new ArrayList<>(results.size());
    for (final DoubletResult result : results) {
        final double score = result.getMaxScore();
        // Filter the singles that would be accepted
        if (result.good1) {
            // Filter the doublets that would be accepted
            if (result.good2) {
                final double tp2 = result.tp2a + result.tp2b;
                final double fp2 = result.fp2a + result.fp2b;
                tp += tp2;
                fp += fp2;
                if (result.tp2a > 0.5) {
                    bestTp += result.tp2a;
                    bestFp += result.fp2a;
                }
                if (result.tp2b > 0.5) {
                    bestTp += result.tp2b;
                    bestFp += result.fp2b;
                }
                // Store this as a doublet bonus
                data.add(new DoubletBonus(score, result.getAvScore(), tp2 - result.tp1, fp2 - result.fp1));
            } else {
                // No doublet fit so this will always be the single fit result
                tp += result.tp1;
                fp += result.fp1;
                if (result.tp1 > 0.5) {
                    bestTp += result.tp1;
                    bestFp += result.fp1;
                }
            }
        }
        // Build statistics
        final int c = result.matchClass;
        // True results, i.e. where there was a choice between single or doublet
        if (result.valid) {
            stats[c].add(score);
        }
        // Of those where the fit was good, summarise the iterations and evaluations
        if (result.good1) {
            stats[3].add(result.iter1);
            stats[4].add(result.eval1);
            // about the iteration increase for singles that are not doublets.
            if (c != 0 && result.good2) {
                stats[5].add(result.iter2);
                stats[6].add(result.eval2);
            }
        }
    }
    // Debug the counts
    // double tpSingle = 0;
    // double fpSingle = 0;
    // double tpDoublet = 0;
    // double fpDoublet = 0;
    // int nSingle = 0, nDoublet = 0;
    // for (DoubletResult result : results) {
    // if (result.good1) {
    // if (result.good2) {
    // tpDoublet += result.tp2a + result.tp2b;
    // fpDoublet += result.fp2a + result.fp2b;
    // nDoublet++;
    // }
    // tpSingle += result.tp1;
    // fpSingle += result.fp1;
    // nSingle++;
    // }
    // }
    // System.out.printf("Single %.1f,%.1f (%d) : Doublet %.1f,%.1f (%d)\n", tpSingle, fpSingle,
    // nSingle, tpDoublet, fpDoublet, nDoublet * 2);
    // Summarise score for true results
    final Percentile p = new Percentile(99);
    for (int c = 0; c < stats.length; c++) {
        final double[] values = stats[c].getValues();
        // Sorting is need for the percentile and the cumulative histogram so do it once
        Arrays.sort(values);
        sb.append(MathUtils.rounded(stats[c].getMean())).append("+/-").append(MathUtils.rounded(stats[c].getStandardDeviation())).append(" (").append(stats[c].getN()).append(") ").append(MathUtils.rounded(p.evaluate(values))).append('\t');
        if (settings.showHistograms && settings.displayHistograms[c + Settings.NAMES.length]) {
            showCumulativeHistogram(values, Settings.NAMES2[c]);
        }
    }
    sb.append(Settings.MATCHING_METHODS[settings.matchingMethod]).append('\t');
    // Plot a graph of the additional results we would fit at all score thresholds.
    // This assumes we just pick the the doublet if we fit it (NO FILTERING at all!)
    // Initialise the score for residuals 0
    // Store this as it serves as a baseline for the filtering analysis
    final ResidualsScore residualsScoreMax = computeScores(data, tp, fp, numberOfMolecules, true);
    final ResidualsScore residualsScoreAv = computeScores(data, tp, fp, numberOfMolecules, false);
    residualsScore = (settings.useMaxResiduals) ? residualsScoreMax : residualsScoreAv;
    if (settings.showJaccardPlot) {
        plotJaccard(residualsScore, null);
    }
    final String bestJaccard = MathUtils.rounded(bestTp / (bestFp + numberOfMolecules)) + '\t';
    final String analysisPrefix2 = analysisPrefix + bestJaccard;
    sb.append(bestJaccard);
    addJaccardScores(sb);
    sb.append('\t').append(TextUtils.nanosToString(runTime));
    createSummaryTable().append(sb.toString());
    // Store results in memory for later analysis
    referenceResults.set(new ReferenceResults(results, residualsScoreMax, residualsScoreAv, numberOfMolecules, analysisPrefix2));
}

Example 2 with CalibrationReader

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

the class PeakFit method checkCameraModel.

/**
 * Check the camera model covers the region of the source.
 *
 * @param fitConfig the fit config
 * @param sourceBounds the source bounds of the input image
 * @param cropBounds the crop bounds (relative to the input image)
 * @param initialise the initialise flag
 * @return true, if successful
 * @throws IllegalStateException if no camera model exists for the camera type
 */
private static boolean checkCameraModel(FitConfiguration fitConfig, Rectangle sourceBounds, Rectangle cropBounds, boolean initialise) {
    final CalibrationReader calibration = fitConfig.getCalibrationReader();
    if (calibration.isScmos() && sourceBounds != null) {
        CameraModel cameraModel = fitConfig.getCameraModel();
        // The camera model origin must be reset to be relative to the source bounds origin
        cameraModel = cropCameraModel(cameraModel, sourceBounds, cropBounds, true);
        if (cameraModel == null) {
            return false;
        }
        if (initialise && cameraModel instanceof PerPixelCameraModel) {
            ((PerPixelCameraModel) cameraModel).initialise();
        }
        fitConfig.setCameraModel(cameraModel);
    }
    return true;
}

Example 3 with CalibrationReader

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

the class Fire method initialise.

/**
 * Initialise this instance with localisation results for the FIRE computation.
 *
 * @param results the results
 * @param results2 the second set of results (can be null)
 */
public void initialise(MemoryPeakResults results, MemoryPeakResults results2) {
    this.results = verify(results);
    this.results2 = verify(results2);
    if (this.results == null) {
        return;
    }
    nmPerUnit = 1;
    DistanceUnit unit = null;
    units = "unknown";
    final CalibrationReader cal = results.getCalibrationReader();
    if (cal != null) {
        try {
            nmPerUnit = cal.getDistanceConverter(DistanceUnit.NM).convert(1);
            units = UnitHelper.getShortName(DistanceUnit.NM);
            unit = DistanceUnit.NM;
        } catch (final ConversionException ex) {
            IJ.log(pluginTitle + " Warning: Ignoring invalid distance calibration for primary results");
        }
    } else {
        IJ.log(pluginTitle + " Warning: No calibration exists for primary results");
    }
    // Calibration must match between datasets
    if (this.results2 != null) {
        CalibrationReader cal2 = results.getCalibrationReader();
        if (unit == null) {
            if (cal2 != null) {
                IJ.log(pluginTitle + " Warning: Ignoring calibration for secondary results since no calibration" + " exists for primary results");
            }
        } else {
            // The calibration must match
            try {
                // Try to create a converter and check it is the same conversion
                if (cal2 != null && cal2.getDistanceConverter(DistanceUnit.NM).convert(1) != nmPerUnit) {
                    // Set to null to mark invalid
                    cal2 = null;
                }
            } catch (final ConversionException ex) {
                // Set to null to mark invalid
                cal2 = null;
            } finally {
                if (cal2 == null) {
                    this.results = null;
                    IJ.error(pluginTitle, "Error: Calibration between the two input datasets does not match");
                }
            }
        }
    }
    // Use the float data bounds. This prevents problems if the data is far from the origin.
    dataBounds = results.getDataBounds(null);
    if (this.results2 != null) {
        final Rectangle2D dataBounds2 = results.getDataBounds(null);
        dataBounds = dataBounds.createUnion(dataBounds2);
    }
}

Example 4 with CalibrationReader

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

the class TsfPeakResultsWriter method createSpotList.

private SpotList createSpotList() {
    final SpotList.Builder builder = SpotList.newBuilder();
    builder.setApplicationId(APPLICATION_ID);
    builder.setNrSpots(size);
    // Add the standard details the TSF supports. We use extensions to add GDSC SMLM data.
    if (!TextUtils.isNullOrEmpty(getName())) {
        builder.setName(getName());
    }
    if (getSource() != null) {
        builder.setNrPixelsX(getSource().width);
        builder.setNrPixelsY(getSource().height);
        builder.setNrFrames(getSource().frames);
        builder.setSource(singleLine(getSource().toXml()));
    }
    if (getBounds() != null) {
        final ROI.Builder roiBuilder = builder.getRoiBuilder();
        roiBuilder.setX(getBounds().x);
        roiBuilder.setY(getBounds().y);
        roiBuilder.setXWidth(getBounds().width);
        roiBuilder.setYWidth(getBounds().height);
        builder.setRoi(roiBuilder.build());
    }
    if (hasCalibration()) {
        final CalibrationReader cr = getCalibrationReader();
        if (cr.hasNmPerPixel()) {
            builder.setPixelSize((float) cr.getNmPerPixel());
        }
        if (cr.hasExposureTime()) {
            builder.setExposureTime(cr.getExposureTime());
        }
        if (cr.hasReadNoise()) {
            builder.setReadNoise(cr.getReadNoise());
        }
        if (cr.hasBias()) {
            builder.setBias(cr.getBias());
        }
        if (cr.hasCameraType()) {
            builder.setCameraType(cameraTypeMap[cr.getCameraType().ordinal()]);
        }
        if (cr.hasDistanceUnit()) {
            builder.setLocationUnits(locationUnitsMap[cr.getDistanceUnit().ordinal()]);
        }
        if (cr.hasIntensityUnit()) {
            builder.setIntensityUnits(intensityUnitsMap[cr.getIntensityUnit().ordinal()]);
        }
        if (cr.hasAngleUnit()) {
            builder.setThetaUnits(thetaUnitsMap[cr.getAngleUnit().ordinal()]);
        }
        // We can use some logic here to get the QE
        if (cr.hasCountPerPhoton()) {
            builder.setGain(cr.getCountPerPhoton());
            final double qe = (cr.hasQuantumEfficiency()) ? cr.getQuantumEfficiency() : 1;
            // e-/photon / count/photon => e-/count
            final double ecf = qe / cr.getCountPerPhoton();
            builder.addEcf(ecf);
            builder.addQe(qe);
        }
    }
    if (!TextUtils.isNullOrEmpty(getConfiguration())) {
        builder.setConfiguration(singleLine(getConfiguration()));
    }
    if (getPsf() != null) {
        try {
            final Printer printer = JsonFormat.printer().omittingInsignificantWhitespace();
            builder.setPSF(printer.print(getPsf()));
        } catch (final InvalidProtocolBufferException ex) {
            // This shouldn't happen so throw it
            throw new NotImplementedException("Unable to serialise the PSF settings", ex);
        }
    }
    // Have a property so the boxSize can be set
    if (boxSize > 0) {
        builder.setBoxSize(boxSize);
    }
    builder.setFitMode(fitMode);
    final FluorophoreType.Builder typeBuilder = FluorophoreType.newBuilder();
    typeBuilder.setId(1);
    typeBuilder.setDescription("Default fluorophore");
    typeBuilder.setIsFiducial(false);
    builder.addFluorophoreTypes(typeBuilder.build());
    return builder.build();
}

Example 5 with CalibrationReader

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

the class MalkFilePeakResults method getHeaderComments.

@Override
protected String[] getHeaderComments() {
    final String[] comments = new String[3];
    int count = 0;
    if (hasCalibration()) {
        final CalibrationReader cr = getCalibrationReader();
        if (cr.hasNmPerPixel()) {
            comments[count++] = String.format("Pixel pitch %s (nm)", MathUtils.rounded(cr.getNmPerPixel()));
        }
        if (cr.hasCountPerPhoton()) {
            comments[count++] = String.format("Gain %s (Count/photon)", MathUtils.rounded(cr.getCountPerPhoton()));
        }
        if (cr.hasExposureTime()) {
            comments[count++] = String.format("Exposure time %s (seconds)", MathUtils.rounded(cr.getExposureTime() * 1e-3));
        }
    }
    return Arrays.copyOf(comments, count);
}