use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class PcPalmFitting method fitRandomModel.
/**
* Fits the correlation curve with r>0 to the random model using the estimated density and
* precision. Parameters must be fit within a tolerance of the starting values.
*
* @param gr the correlation curve
* @param sigmaS The estimated precision
* @param proteinDensity The estimate protein density
* @param resultColour the result colour
* @return The fitted parameters [precision, density]
*/
@Nullable
private double[] fitRandomModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
final RandomModelFunction function = new RandomModelFunction();
randomModel = function;
ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", randomModel.getName(), sigmaS, proteinDensity * 1e6);
randomModel.setLogging(true);
for (int i = offset; i < gr[0].length; i++) {
// error)
if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
randomModel.addPoint(gr[0][i], gr[1][i]);
}
}
final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
Optimum optimum;
try {
// @formatter:off
final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(new double[] { sigmaS, proteinDensity }).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
// @formatter:on
optimum = optimizer.optimize(problem);
} catch (final TooManyIterationsException ex) {
ImageJUtils.log("Failed to fit %s: Too many iterations (%s)", randomModel.getName(), ex.getMessage());
return null;
} catch (final ConvergenceException ex) {
ImageJUtils.log("Failed to fit %s: %s", randomModel.getName(), ex.getMessage());
return null;
}
randomModel.setLogging(false);
final double[] parameters = optimum.getPoint().toArray();
// Ensure the width is positive
parameters[0] = Math.abs(parameters[0]);
final double ss = optimum.getResiduals().dotProduct(optimum.getResiduals());
final double totalSumSquares = MathUtils.getTotalSumOfSquares(randomModel.getY());
randomModelAdjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, randomModel.size(), parameters.length);
final double fitSigmaS = parameters[0];
final double fitProteinDensity = parameters[1];
// Check the fitted parameters are within tolerance of the initial estimates
final double e1 = parameterDrift(sigmaS, fitSigmaS);
final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
ImageJUtils.log(" %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", randomModel.getName(), ss, randomModelAdjustedR2, optimum.getEvaluations());
ImageJUtils.log(" %s parameters:", randomModel.getName());
ImageJUtils.log(" Average precision = %s nm (%s%%)", MathUtils.rounded(fitSigmaS, 4), MathUtils.rounded(e1, 4));
ImageJUtils.log(" Average protein density = %s um^-2 (%s%%)", MathUtils.rounded(fitProteinDensity * 1e6, 4), MathUtils.rounded(e2, 4));
valid1 = true;
if (settings.fittingTolerance > 0 && (Math.abs(e1) > settings.fittingTolerance || Math.abs(e2) > settings.fittingTolerance)) {
ImageJUtils.log(" Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%)," + " average protein density = %g um^-2 (%s%%)", randomModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
valid1 = false;
}
if (valid1) {
// ---------
// TODO - My data does not comply with this criteria.
// This could be due to the PC-PALM Molecule code limiting the nmPerPixel to fit the images in
// memory thus removing correlations at small r.
// It could also be due to the nature of the random simulations being 3D not 2D membranes
// as per the PC-PALM paper.
// ---------
// Evaluate g(r)protein where:
// g(r)peaks = g(r)protein + g(r)stoch
// g(r)peaks ~ 1 + g(r)stoch
// Verify g(r)protein should be <1.5 for all r>0
final double[] grStoch = randomModel.value(parameters);
final double[] grPeaks = randomModel.getY();
final double[] radius = randomModel.getX();
for (int i = 0; i < grPeaks.length; i++) {
// Only evaluate above the fitted average precision
if (radius[i] < fitSigmaS) {
continue;
}
// Note the RandomModelFunction evaluates g(r)stoch + 1
final double grProtein = grPeaks[i] - (grStoch[i] - 1);
if (grProtein > settings.grProteinThreshold) {
// Failed fit
ImageJUtils.log(" Failed to fit %s: g(r)protein %s > %s @ r=%s", randomModel.getName(), MathUtils.rounded(grProtein, 4), MathUtils.rounded(settings.grProteinThreshold, 4), MathUtils.rounded(radius[i], 4));
valid1 = false;
}
}
}
addResult(randomModel.getName(), resultColour, valid1, fitSigmaS, fitProteinDensity, 0, 0, 0, 0, randomModelAdjustedR2);
return parameters;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class PcPalmMolecules method plot.
@Nullable
private static double[][] plot(DoubleData stats, String label, boolean integerBins) {
if (integerBins) {
// The histogram is not need for the return statement
new HistogramPlotBuilder(TITLE, stats, label).setMinBinWidth(1).show();
return null;
}
// Show a cumulative histogram so that the bin size is not relevant
final double[][] hist = MathUtils.cumulativeHistogram(stats.values(), false);
// Create the axes
final double[] xValues = hist[0];
final double[] yValues = hist[1];
// Plot
final String title = TITLE + " " + label;
final Plot plot = new Plot(title, label, "Frequency");
plot.addPoints(xValues, yValues, Plot.LINE);
ImageJUtils.display(title, plot);
return hist;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class PeakResultsReader method createNStormResult.
/**
* Creates the NStorm result.
*
* @param line the line
* @param readPhotons Set to {@code true} if there is a Photons field
* @return the peak result
*/
@Nullable
private static PeakResult createNStormResult(String line, boolean readPhotons) {
// Ywc
try (Scanner scanner = new Scanner(line)) {
scanner.useDelimiter(tabPattern);
scanner.useLocale(Locale.US);
// channelName
scanner.next();
// x
scanner.nextFloat();
// y
scanner.nextFloat();
final float xc = scanner.nextFloat();
final float yc = scanner.nextFloat();
final float height = scanner.nextFloat();
final float area = scanner.nextFloat();
final float width = scanner.nextFloat();
// phi
scanner.nextFloat();
final float ax = scanner.nextFloat();
final float bg = scanner.nextFloat();
// intensity
scanner.nextFloat();
final int frame = scanner.nextInt();
final int length = scanner.nextInt();
double photons = 0;
if (readPhotons) {
// These are not needed
// Link
scanner.next();
// Valid
scanner.next();
// Z
scanner.next();
// Zc
scanner.next();
photons = scanner.nextDouble();
}
// The coordinates are in nm
// The values are in ADUs. The area value is the signal.
// The following relationship holds when length == 1:
// Area = Height * 2 * pi * (Width / (pixel_pitch*2) )^2
// => Pixel_pitch = 0.5 * Width / sqrt(Area / (Height * 2 * pi))
final float[] params = new float[SIZE_TWO_AXIS];
params[PeakResult.BACKGROUND] = bg;
params[PeakResult.INTENSITY] = area;
params[PeakResult.X] = xc;
params[PeakResult.Y] = yc;
// Convert width (2*SD) to SD
final float sd = width / 2f;
// Convert to separate XY widths using the axial ratio
if (ax == 1) {
params[INDEX_SX] = sd;
params[INDEX_SY] = sd;
} else {
// Ensure the axial ratio is long/short
final double a = Math.sqrt((ax < 1) ? 1.0 / ax : ax);
params[INDEX_SX] = (float) (sd * a);
params[INDEX_SY] = (float) (sd / a);
}
// Store the photons in the error value
return new ExtendedPeakResult(frame, (int) xc, (int) yc, height, photons, 0.0f, 0, params, null, frame + length - 1, 0);
} catch (final NoSuchElementException ignored) {
// Ignore
}
return null;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class BenchmarkFilterAnalysis method score.
@Nullable
@Override
public SearchResult<FilterScore>[] score(double[][] points) {
gaIteration++;
SimpleFilterScore max = filterScoreOptimum;
final FilterScoreResult[] scoreResults = scoreFilters(setStrength(new FilterSet(searchSpaceToFilters(points))), false);
if (scoreResults == null) {
return null;
}
@SuppressWarnings("unchecked") final SearchResult<FilterScore>[] scores = new SearchResult[scoreResults.length];
for (int index = 0; index < scoreResults.length; index++) {
final FilterScoreResult scoreResult = scoreResults[index];
final SimpleFilterScore result = new SimpleFilterScore(scoreResult, true, scoreResult.criteria >= minCriteria);
if (result.compareTo(max) < 0) {
max = result;
}
scores[index] = new SearchResult<>(result.result.filter.getParameters(), result);
}
filterScoreOptimum = max;
// Add the best filter to the table
// This filter may not have been part of the scored subset so use the entire results set for
// reporting
final DirectFilter filter = max.result.filter;
final FractionClassificationResult r = scoreFilter(filter, DEFAULT_MINIMUM_FILTER, gaResultsList, coordinateStore);
final StringBuilder text = createResult(filter, r);
add(text, gaIteration);
gaWindow.accept(text.toString());
return scores;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.
the class BenchmarkFilterAnalysis method readFilterSets.
@Nullable
@SuppressWarnings("unchecked")
private List<FilterSet> readFilterSets() {
if (extraOptions) {
final MultiPathFilter multiFilter = BenchmarkSpotFit.getMultiFilter();
if (multiFilter != null) {
final IDirectFilter f = multiFilter.getFilter();
if (f instanceof DirectFilter) {
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage("Use an identical filter to " + BenchmarkSpotFit.TITLE);
gd.enableYesNoCancel();
gd.hideCancelButton();
gd.showDialog();
if (gd.wasOKed()) {
final List<FilterSet> filterSets = new ArrayList<>(1);
final List<Filter> filters = new ArrayList<>(1);
filters.add((DirectFilter) f);
final FilterSet filterSet = new FilterSet(filters);
filterSets.add(filterSet);
resetParametersFromFitting();
createResultsPrefix2();
return filterSets;
}
}
}
}
GUIFilterSettings filterSettings = SettingsManager.readGuiFilterSettings(0);
final String filename = ImageJUtils.getFilename("Filter_File", filterSettings.getFilterSetFilename());
if (filename != null) {
IJ.showStatus("Reading filters ...");
filterSettings = filterSettings.toBuilder().setFilterSetFilename(filename).build();
// Allow the filters to be cached
final Triple<String, Long, List<FilterSet>> filterCache = LAST_FILTER_LIST.get();
if (isSameFile(filename, filterCache)) {
final GenericDialog gd = new GenericDialog(TITLE);
gd.hideCancelButton();
gd.addMessage("The same filter file was selected.");
gd.addCheckbox("Re-use_filters", settings.reUseFilters);
gd.showDialog();
if (!gd.wasCanceled()) {
settings.reUseFilters = gd.getNextBoolean();
if (settings.reUseFilters) {
SettingsManager.writeSettings(filterSettings);
return filterCache.getRight();
}
}
}
final Path path = Paths.get(filename);
try (BufferedReader input = new BufferedReader(new UnicodeReader(Files.newInputStream(path), null))) {
// Use the instance so we can catch the exception
final Object o = FilterXStreamUtils.getXStreamInstance().fromXML(input);
if (!(o instanceof List<?>)) {
IJ.log("No filter sets defined in the specified file: " + filename);
return null;
}
SettingsManager.writeSettings(filterSettings);
List<FilterSet> filterSets = (List<FilterSet>) o;
if (containsStandardFilters(filterSets)) {
IJ.log("Filter sets must contain 'Direct' filters");
return null;
}
// Check they are not empty lists
final List<FilterSet> filterSets2 = new LinkedList<>();
for (final FilterSet filterSet : filterSets) {
if (filterSet.size() != 0) {
filterSets2.add(filterSet);
} else {
IJ.log("Filter set empty: " + filterSet.getName());
}
}
if (filterSets2.isEmpty()) {
IJ.log("All Filter sets are empty");
return null;
}
// Maintain the same list type
filterSets.clear();
filterSets.addAll(filterSets2);
// Option to enumerate filters
filterSets = expandFilters(filterSets);
// Save for re-use
LAST_FILTER_LIST.set(Triple.of(filename, getLastModified(path), filterSets));
return filterSets;
} catch (final Exception ex) {
IJ.log("Unable to load the filter sets from file: " + ex.getMessage());
} finally {
IJ.showStatus("");
}
}
return null;
}
Aggregations