Java tutorial
package gdsc.smlm.ij.plugins; /*----------------------------------------------------------------------------- * GDSC SMLM Software * * Copyright (C) 2013 Alex Herbert * Genome Damage and Stability Centre * University of Sussex, UK * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 3 of the License, or * (at your option) any later version. *---------------------------------------------------------------------------*/ import gdsc.smlm.fitting.FitResult; import gdsc.smlm.fitting.FitStatus; import gdsc.smlm.function.gaussian.Gaussian2DFunction; import gdsc.smlm.ij.plugins.BenchmarkSpotFit.FilterCandidates; import gdsc.smlm.ij.settings.FilterSettings; import gdsc.smlm.ij.settings.GlobalSettings; import gdsc.smlm.ij.settings.SettingsManager; import gdsc.smlm.ij.utils.Utils; import gdsc.smlm.results.Calibration; import gdsc.smlm.results.MemoryPeakResults; import gdsc.smlm.results.PeakResult; import gdsc.smlm.results.filter.Filter; import gdsc.smlm.results.filter.FilterSet; import gdsc.smlm.results.filter.XStreamWrapper; import gdsc.smlm.results.match.FractionClassificationResult; import gdsc.smlm.utils.Maths; import gdsc.smlm.utils.Sort; import gdsc.smlm.utils.StoredDataStatistics; import gdsc.smlm.utils.UnicodeReader; import gdsc.smlm.utils.XmlUtils; import ij.IJ; import ij.gui.GenericDialog; import ij.gui.Plot2; import ij.gui.PlotWindow; import ij.plugin.PlugIn; import ij.plugin.WindowOrganiser; import ij.text.TextWindow; import java.awt.Color; import java.awt.Point; import java.io.BufferedReader; import java.io.FileInputStream; import java.io.FileOutputStream; import java.io.IOException; import java.io.OutputStreamWriter; import java.util.ArrayList; import java.util.Arrays; import java.util.Collections; import java.util.HashMap; import java.util.LinkedList; import java.util.List; import java.util.Map.Entry; import org.apache.commons.math3.analysis.interpolation.LoessInterpolator; import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; /** * Run different filtering methods on a set of benchmark fitting results outputting performance statistics on the * success of the filter. The fitting results are generated by the BenchmarkSpotFit plugin. * <p> * Filtering is done using e.g. SNR threshold, Precision thresholds, etc. The statistics reported are shown in a table, * e.g. precision, Jaccard, F-score. */ public class BenchmarkFilterAnalysis implements PlugIn { private static final String TITLE = "Filter Analysis"; private static TextWindow resultsWindow = null; private static TextWindow summaryWindow = null; private static TextWindow sensitivityWindow = null; private static int failCount = 1; private static int failCountRange = 0; private static boolean rerankBySignal = false; private static boolean showResultsTable = false; private static boolean showSummaryTable = true; private static boolean clearTables = false; private static String filterFilename = ""; private static int summaryTopN = 0; private static int plotTopN = 0; private static boolean saveBestFilter = false; private ArrayList<NamedPlot> plots; private static boolean calculateSensitivity = false; private static double delta = 0.1; private static int criteriaIndex; private static double criteriaLimit = 0.9; private double minCriteria = 0; private boolean invertCriteria = false; private static int scoreIndex; private boolean invertScore = false; private static double upperMatchDistance = 100; private static double partialMatchDistance = 100; private static boolean depthRecallAnalysis = true; private static String resultsTitle; private String resultsPrefix, resultsPrefix2; private static String resultsPrefix3; private HashMap<String, FilterScore> bestFilter; private LinkedList<String> bestFilterOrder; private static boolean reUseFilters = true; private static String oldFilename = ""; private static List<FilterSet> filterList = null; private static int lastId = 0; private static boolean lastRank = false; private static double lastUpperMatchDistance = -1, lastPartialMatchDistance = -1; private static List<MemoryPeakResults> resultsList = null; private static StoredDataStatistics depthStats, depthFitStats; private boolean isHeadless; private static String[] COLUMNS = { "nP", "TP", "FP", "TN", "FN", "TPR", "TNR", "PPV", "NPV", "FPR", "FNR", "FDR", "ACC", "MCC", "Informedness", "Markedness", "Recall", "Precision", "F1", "Jaccard" }; private static boolean[] showColumns; static { showColumns = new boolean[COLUMNS.length]; Arrays.fill(showColumns, true); showColumns[0] = false; // nP criteriaIndex = COLUMNS.length - 3; scoreIndex = COLUMNS.length - 1; } private CreateData.SimulationParameters simulationParameters; public BenchmarkFilterAnalysis() { isHeadless = java.awt.GraphicsEnvironment.isHeadless(); } /* * (non-Javadoc) * * @see ij.plugin.PlugIn#run(java.lang.String) */ public void run(String arg) { simulationParameters = CreateData.simulationParameters; if (simulationParameters == null) { IJ.error(TITLE, "No benchmark spot parameters in memory"); return; } if (BenchmarkSpotFit.fitResults == null) { IJ.error(TITLE, "No benchmark fitting results in memory"); return; } if (BenchmarkSpotFit.lastId != simulationParameters.id) { IJ.error(TITLE, "Update the benchmark spot fitting for the latest simulation"); return; } if (BenchmarkSpotFit.lastFilterId != BenchmarkSpotFilter.filterResultsId) { IJ.error(TITLE, "Update the benchmark spot fitting for the latest filter"); return; } resultsList = readResults(); if (resultsList.isEmpty()) { IJ.error(TITLE, "No results could be loaded"); return; } if (!showDialog(resultsList)) return; // Load filters from file List<FilterSet> filterSets = readFilterSets(); if (filterSets == null || filterSets.isEmpty()) { IJ.error(TITLE, "No filters specified"); return; } long time = analyse(resultsList, filterSets); IJ.showStatus("Finished : " + Utils.timeToString(time)); } @SuppressWarnings("unchecked") private List<FilterSet> readFilterSets() { GlobalSettings gs = SettingsManager.loadSettings(); FilterSettings filterSettings = gs.getFilterSettings(); String filename = Utils.getFilename("Filter_File", filterSettings.filterSetFilename); if (filename != null) { IJ.showStatus("Reading filters ..."); filterSettings.filterSetFilename = filename; // Allow the filters to be cached if (filename.equals(oldFilename) && filterList != null) { GenericDialog gd = new GenericDialog(TITLE); gd.hideCancelButton(); gd.addMessage("The same filter file was selected."); gd.addCheckbox("Re-use_filters", reUseFilters); gd.showDialog(); if (!gd.wasCanceled()) { if ((reUseFilters = gd.getNextBoolean())) return filterList; } } BufferedReader input = null; oldFilename = ""; try { FileInputStream fis = new FileInputStream(filterSettings.filterSetFilename); input = new BufferedReader(new UnicodeReader(fis, null)); Object o = XStreamWrapper.getInstance().fromXML(input); if (o != null && o instanceof List<?>) { SettingsManager.saveSettings(gs); filterList = (List<FilterSet>) o; oldFilename = filename; return filterList; } IJ.log("No filter sets defined in the specified file: " + filterSettings.filterSetFilename); } catch (Exception e) { IJ.log("Unable to load the filter sets from file: " + e.getMessage()); } finally { if (input != null) { try { input.close(); } catch (IOException e) { // Ignore } } IJ.showStatus(""); } } return null; } private List<MemoryPeakResults> readResults() { if (resultsList == null || lastId != BenchmarkSpotFit.fitResultsId || lastRank != rerankBySignal || lastPartialMatchDistance != partialMatchDistance || lastUpperMatchDistance != upperMatchDistance) { lastId = BenchmarkSpotFit.fitResultsId; lastRank = rerankBySignal; lastUpperMatchDistance = upperMatchDistance; lastPartialMatchDistance = partialMatchDistance; resultsList = new LinkedList<MemoryPeakResults>(); depthStats = new StoredDataStatistics(); depthFitStats = new StoredDataStatistics(); final double fuzzyMax = BenchmarkSpotFit.distanceInPixels * upperMatchDistance / 100.0; final double fuzzyMin = BenchmarkSpotFit.distanceInPixels * partialMatchDistance / 100.0; final double fuzzyRange = fuzzyMax - fuzzyMin; resultsPrefix3 = "\t" + Utils.rounded(fuzzyMin * simulationParameters.a) + "\t" + Utils.rounded(fuzzyMax * simulationParameters.a); // Normalise the matches so the fuzzy weighting sums to the original true positive count //int tp = 0; //double sum = 0; MemoryPeakResults r = new MemoryPeakResults(); Calibration cal = new Calibration(simulationParameters.a, simulationParameters.gain, 100); cal.bias = simulationParameters.bias; cal.emCCD = simulationParameters.emCCD; cal.readNoise = simulationParameters.readNoise; r.setCalibration(cal); // Set the configuration used for fitting r.setConfiguration(XmlUtils.toXML(BenchmarkSpotFit.fitConfig)); for (Entry<Integer, FilterCandidates> entry : BenchmarkSpotFit.fitResults.entrySet()) { final int peak = entry.getKey().intValue(); final FilterCandidates result = entry.getValue(); depthStats.add(result.zPosition); // Results are in order of candidate ranking. int[] indices = Utils.newArray(result.fitResult.length, 0, 1); if (rerankBySignal) { // Change to signal intensity. double[] score = new double[result.fitResult.length]; for (int i = 0; i < result.fitResult.length; i++) { final FitResult fitResult = result.fitResult[i]; if (fitResult.getStatus() == FitStatus.OK) { score[i] = fitResult.getParameters()[Gaussian2DFunction.SIGNAL]; } } Sort.sort(indices, score); } // Data about the match is stored in small arrays of size equal to the number of matches. // Create an array holding the full index position corresponding to each match. int[] positionIndex = new int[result.dMatch.length]; int count = 0; for (int i = 0; i < result.fitMatch.length; i++) { if (result.fitMatch[i]) positionIndex[count++] = i; } for (int i = 0; i < result.fitResult.length; i++) { final int index = indices[i]; final FitResult fitResult = result.fitResult[index]; if (fitResult.getStatus() == FitStatus.OK) { // Assume we are not fitting doublets and the fit result will have 1 peak final float[] params = Utils.toFloat(fitResult.getParameters()); // To allow the shift filter to function the X/Y position coordinates must be relative final double[] initial = fitResult.getInitialParameters(); params[Gaussian2DFunction.X_POSITION] -= initial[Gaussian2DFunction.X_POSITION]; params[Gaussian2DFunction.Y_POSITION] -= initial[Gaussian2DFunction.Y_POSITION]; // These will be incorrect due to the relative adjustment above... final int origX = (int) params[Gaussian2DFunction.X_POSITION]; final int origY = (int) params[Gaussian2DFunction.Y_POSITION]; // Binary classification uses a score of 1 (match) or 0 (no match) // Fuzzy classification uses a score of 1 (match) or 0 (no match), or >0 <1 (partial match) float score = 0; double depth = 0; if (result.fitMatch[index]) { int position = getPosition(positionIndex, index); double distance = result.dMatch[position]; // Store depth of matches for later analysis depth = result.zMatch[position]; if (distance <= fuzzyMax) { depthFitStats.add(depth); if (distance <= fuzzyMin) { score = 1; } else { // Interpolate from the minimum to the maximum match distance: // Linear //score = (float) ((fuzzyMax - result.d[index]) / fuzzyRange); // Cosine score = (float) (0.5 * (1 + Math.cos(((distance - fuzzyMin) / fuzzyRange) * Math.PI))); } //tp++; //sum += score; } } // Use a custom peak result so that the depth can be stored. r.add(new DepthPeakResult(peak, origX, origY, score, fitResult.getError(), result.noise, params, depth)); } } } if (r.size() > 0) { // Normalise the fuzzy weighting //final double w = tp / sum; //for (PeakResult result : r.getResults()) //{ // if (result.origValue != 0) // result.origValue *= w; //} resultsList.add(r); } } return resultsList; } private int getPosition(int[] positionIndex, int index) { for (int c = 0; c < positionIndex.length; c++) if (positionIndex[c] == index) return c; throw new RuntimeException("Cannot find the position for the match"); } private boolean showDialog(List<MemoryPeakResults> resultsList) { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); int total = 0; int tp = 0; for (MemoryPeakResults r : resultsList) { total += r.size(); for (PeakResult p : r.getResults()) if (p.origValue != 0) tp++; } gd.addMessage(String.format("%d results, %d True-Positives", total, tp)); gd.addSlider("Fail_count", 0, 20, failCount); gd.addSlider("Fail_count_range", 0, 5, failCountRange); gd.addCheckbox("Rank_by_signal", rerankBySignal); gd.addCheckbox("Show_table", showResultsTable); gd.addCheckbox("Show_summary", showSummaryTable); gd.addCheckbox("Clear_tables", clearTables); gd.addSlider("Summary_top_n", 0, 20, summaryTopN); gd.addSlider("Plot_top_n", 0, 20, plotTopN); gd.addCheckbox("Save_best_filter", saveBestFilter); gd.addCheckbox("Calculate_sensitivity", calculateSensitivity); gd.addSlider("Delta", 0.01, 1, delta); gd.addChoice("Criteria", COLUMNS, COLUMNS[criteriaIndex]); gd.addNumericField("Criteria_limit", criteriaLimit, 4); gd.addChoice("Score", COLUMNS, COLUMNS[scoreIndex]); gd.addMessage("Fitting results match distance = " + Utils.rounded(BenchmarkSpotFit.distanceInPixels * simulationParameters.a) + " nm"); gd.addSlider("Upper_match_distance (%)", 0, 100, upperMatchDistance); gd.addSlider("Partial_match_distance (%)", 0, 100, partialMatchDistance); if (!simulationParameters.fixedDepth) gd.addCheckbox("Depth_recall_analysis", depthRecallAnalysis); gd.addStringField("Title", resultsTitle, 20); gd.showDialog(); if (gd.wasCanceled() || !readDialog(gd)) return false; if (showResultsTable || showSummaryTable) { gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); gd.addMessage("Select the results:"); for (int i = 0; i < COLUMNS.length; i++) gd.addCheckbox(COLUMNS[i], showColumns[i]); gd.showDialog(); if (gd.wasCanceled()) return false; for (int i = 0; i < COLUMNS.length; i++) showColumns[i] = gd.getNextBoolean(); } // We may have to read the results again if the ranking option has changed if (lastRank != rerankBySignal || lastPartialMatchDistance != partialMatchDistance || lastUpperMatchDistance != upperMatchDistance) readResults(); return true; } private boolean readDialog(GenericDialog gd) { failCount = (int) Math.abs(gd.getNextNumber()); failCountRange = (int) Math.abs(gd.getNextNumber()); rerankBySignal = gd.getNextBoolean(); showResultsTable = gd.getNextBoolean(); showSummaryTable = gd.getNextBoolean(); clearTables = gd.getNextBoolean(); summaryTopN = (int) Math.abs(gd.getNextNumber()); plotTopN = (int) Math.abs(gd.getNextNumber()); saveBestFilter = gd.getNextBoolean(); calculateSensitivity = gd.getNextBoolean(); delta = gd.getNextNumber(); criteriaIndex = gd.getNextChoiceIndex(); criteriaLimit = gd.getNextNumber(); scoreIndex = gd.getNextChoiceIndex(); upperMatchDistance = Math.abs(gd.getNextNumber()); partialMatchDistance = Math.abs(gd.getNextNumber()); if (!simulationParameters.fixedDepth) depthRecallAnalysis = gd.getNextBoolean(); resultsTitle = gd.getNextString(); resultsPrefix = BenchmarkSpotFit.resultPrefix + "\t" + resultsTitle + "\t"; resultsPrefix2 = "\t" + failCount; if (failCountRange > 0) resultsPrefix2 += "-" + (failCount + failCountRange); // Check there is one output if (!showResultsTable && !showSummaryTable && !calculateSensitivity && plotTopN < 1 && !saveBestFilter) { IJ.error(TITLE, "No output selected"); return false; } // Check arguments try { Parameters.isAboveZero("Delta", delta); Parameters.isBelow("Delta", delta, 1); Parameters.isEqualOrBelow("Partial match distance", partialMatchDistance, upperMatchDistance); } catch (IllegalArgumentException e) { IJ.error(TITLE, e.getMessage()); return false; } invertCriteria = requiresInversion(criteriaIndex); minCriteria = (invertCriteria) ? -criteriaLimit : criteriaLimit; invertScore = requiresInversion(scoreIndex); return !gd.invalidNumber(); } /** * Run different filtering methods on a set of labelled peak results outputting performance statistics on the * success of * the filter to an ImageJ table. * <p> * If the peak result original value is set to 1 it is considered a true peak, 0 for a false peak. Filtering is done * using e.g. SNR threshold, Precision thresholds, etc. The statistics reported are shown in a table, e.g. * precision, Jaccard, F-score. * <p> * For each filter set a plot is shown of the score verses the filter value, thus filters should be provided in * ascending numerical order otherwise they are sorted. * * @param resultsList * @param filterSets */ private long analyse(List<MemoryPeakResults> resultsList, List<FilterSet> filterSets) { createResultsWindow(); plots = new ArrayList<NamedPlot>(plotTopN); bestFilter = new HashMap<String, FilterScore>(); bestFilterOrder = new LinkedList<String>(); IJ.showStatus("Analysing filters ..."); long time = System.currentTimeMillis(); int total = countFilters(filterSets); int count = 0; for (FilterSet filterSet : filterSets) { IJ.showStatus("Analysing " + filterSet.getName() + " ..."); count = run(filterSet, resultsList, count, total); } time = System.currentTimeMillis() - time; IJ.showProgress(1); IJ.showStatus(""); if (bestFilter.isEmpty()) { IJ.log("Warning: No filters pass the criteria"); return time; } List<FilterScore> filters = new ArrayList<FilterScore>(bestFilter.values()); if (showSummaryTable || saveBestFilter) Collections.sort(filters); if (showSummaryTable) { createSummaryWindow(); int n = 0; for (FilterScore fs : filters) { FractionClassificationResult r = scoreFilter(fs.filter, resultsList); String text = createResult(fs.filter, r); if (isHeadless) IJ.log(text); else summaryWindow.append(text); n++; if (summaryTopN > 0 && n >= summaryTopN) break; } // Add a spacer to the summary table if we have multiple results if (n > 1) { if (isHeadless) IJ.log(""); else summaryWindow.append(""); } } if (saveBestFilter) saveFilter(filters.get(0).filter); showPlots(); calculateSensitivity(resultsList); depthAnalysis(filters.get(0).filter); return time; } private int countFilters(List<FilterSet> filterSets) { int count = 0; for (FilterSet filterSet : filterSets) count += filterSet.size(); return count; } private void showPlots() { if (plots.isEmpty()) return; // Display the top N plots int[] list = new int[plots.size()]; int i = 0; for (NamedPlot p : plots) { Plot2 plot = new Plot2(p.name, p.xAxisName, COLUMNS[scoreIndex], p.xValues, p.yValues); plot.setLimits(p.xValues[0], p.xValues[p.xValues.length - 1], 0, 1); plot.setColor(Color.RED); plot.draw(); plot.setColor(Color.BLUE); plot.addPoints(p.xValues, p.yValues, Plot2.CROSS); PlotWindow plotWindow = Utils.display(p.name, plot); list[i++] = plotWindow.getImagePlus().getID(); } new WindowOrganiser().tileWindows(list); } private void calculateSensitivity(List<MemoryPeakResults> resultsList) { if (!calculateSensitivity) return; if (!bestFilter.isEmpty()) { IJ.showStatus("Calculating sensitivity ..."); createSensitivityWindow(); int currentIndex = 0; for (String type : bestFilterOrder) { IJ.showProgress(currentIndex++, bestFilter.size()); Filter filter = bestFilter.get(type).filter; FractionClassificationResult s = scoreFilter(filter, resultsList); String message = type + "\t\t\t" + Utils.rounded(s.getJaccard(), 4) + "\t\t" + Utils.rounded(s.getPrecision(), 4) + "\t\t" + Utils.rounded(s.getRecall(), 4); if (isHeadless) { IJ.log(message); } else { sensitivityWindow.append(message); } // List all the parameters that can be adjusted. final int parameters = filter.getNumberOfParameters(); for (int index = 0; index < parameters; index++) { // For each parameter compute as upward + downward delta and get the average gradient Filter higher = filter.adjustParameter(index, delta); Filter lower = filter.adjustParameter(index, -delta); FractionClassificationResult sHigher = scoreFilter(higher, resultsList); FractionClassificationResult sLower = scoreFilter(lower, resultsList); StringBuilder sb = new StringBuilder(); sb.append("\t").append(filter.getParameterName(index)).append("\t"); sb.append(Utils.rounded(filter.getParameterValue(index), 4)).append("\t"); double dx1 = higher.getParameterValue(index) - filter.getParameterValue(index); double dx2 = filter.getParameterValue(index) - lower.getParameterValue(index); addSensitivityScore(sb, s.getJaccard(), sHigher.getJaccard(), sLower.getJaccard(), dx1, dx2); addSensitivityScore(sb, s.getPrecision(), sHigher.getPrecision(), sLower.getPrecision(), dx1, dx2); addSensitivityScore(sb, s.getRecall(), sHigher.getRecall(), sLower.getRecall(), dx1, dx2); if (isHeadless) { IJ.log(sb.toString()); } else { sensitivityWindow.append(sb.toString()); } } } String message = "-=-=-=-"; if (isHeadless) { IJ.log(message); } else { sensitivityWindow.append(message); } IJ.showProgress(1); IJ.showStatus(""); } } private void addSensitivityScore(StringBuilder sb, double s, double s1, double s2, double dx1, double dx2) { // Use absolute in case this is not a local maximum. We are mainly interested in how // flat the curve is at this point in relation to parameter changes. double abs1 = Math.abs(s - s1); double abs2 = Math.abs(s - s2); double dydx1 = (abs1) / dx1; double dydx2 = (abs2) / dx2; double relativeSensitivity = (abs1 + abs2) * 0.5; double sensitivity = (dydx1 + dydx2) * 0.5; sb.append(Utils.rounded(relativeSensitivity, 4)).append("\t"); sb.append(Utils.rounded(sensitivity, 4)).append("\t"); } private void createResultsWindow() { if (!showResultsTable) return; if (isHeadless) { IJ.log(createResultsHeader()); } else { if (resultsWindow == null || !resultsWindow.isShowing()) { String header = createResultsHeader(); resultsWindow = new TextWindow(TITLE + " Results", header, "", 900, 300); } if (clearTables) resultsWindow.getTextPanel().clear(); } } private void createSummaryWindow() { if (!showSummaryTable) return; if (isHeadless) { IJ.log(createResultsHeader()); } else { if (summaryWindow == null || !summaryWindow.isShowing()) { String header = createResultsHeader(); summaryWindow = new TextWindow(TITLE + " Summary", header, "", 900, 300); } if (clearTables) summaryWindow.getTextPanel().clear(); } } private String createResultsHeader() { StringBuilder sb = new StringBuilder(BenchmarkSpotFit.tablePrefix); sb.append("\tTitle\tName\tFail\tLower (nm)\tUpper (nm)"); for (int i = 0; i < COLUMNS.length; i++) if (showColumns[i]) sb.append("\t").append(COLUMNS[i]); // Always show the results compared to the original simulation sb.append("\toRecall"); sb.append("\toPrecision"); sb.append("\toF1"); sb.append("\toJaccard"); return sb.toString(); } private void createSensitivityWindow() { if (isHeadless) { IJ.log(createSensitivityHeader()); } else { if (sensitivityWindow == null || !sensitivityWindow.isShowing()) { String header = createSensitivityHeader(); sensitivityWindow = new TextWindow(TITLE + " Sensitivity", header, "", 900, 300); } } } private String createSensitivityHeader() { StringBuilder sb = new StringBuilder(); sb.append("Filter\t"); sb.append("Param\t"); sb.append("Value\t"); sb.append("J Sensitivity (delta)\t"); sb.append("J Sensitivity (unit)\t"); sb.append("P Sensitivity (delta)\t"); sb.append("P Sensitivity (unit)\t"); sb.append("R Sensitivity (delta)\t"); sb.append("R Sensitivity (unit)\t"); return sb.toString(); } private int run(FilterSet filterSet, List<MemoryPeakResults> resultsList, int count, final int total) { double[] xValues = (isHeadless) ? null : new double[filterSet.size()]; double[] yValues = (isHeadless) ? null : new double[filterSet.size()]; int i = 0; filterSet.sort(); // Track if all the filters are the same type. If so then we can calculate the sensitivity of each parameter. String type = null; boolean allSameType = true; Filter maxFilter = null, criteriaFilter = null; double maxScore = -1; double maxCriteria = 0; for (Filter filter : filterSet.getFilters()) { if (count++ % 16 == 0) IJ.showProgress(count, total); if (type == null) type = filter.getType(); else if (!type.equals(filter.getType())) allSameType = false; final double[] result = run(filter, resultsList); final double score = result[0]; final double criteria = result[1]; // Check if the criteria are achieved if (criteria >= minCriteria) { // Check if the score is better if (filter == null || maxScore < score) { maxScore = score; maxFilter = filter; } } else if (maxCriteria < criteria) { maxCriteria = criteria; criteriaFilter = filter; } if (!isHeadless) { xValues[i] = filter.getNumericalValue(); yValues[i++] = score; } } // We may have no filters that pass the criteria if (maxFilter == null) { if (allSameType) { if (criteriaFilter != null) Utils.log("Warning: Filter does not pass the criteria: %s : Best = %s using %s", type, Utils.rounded((invertCriteria) ? -maxCriteria : maxCriteria), criteriaFilter.getName()); else Utils.log("Warning: Filter does not pass the criteria: %s", type); } return count; } if (allSameType) { if (bestFilter.containsKey(type)) { FilterScore filterScore = bestFilter.get(type); if (filterScore.score < maxScore) filterScore.update(maxFilter, maxScore); } else { bestFilter.put(type, new FilterScore(maxFilter, maxScore)); bestFilterOrder.add(type); } } // Add spacer at end of each result set if (isHeadless) { if (showResultsTable && filterSet.size() > 1) IJ.log(""); } else { if (showResultsTable && filterSet.size() > 1) resultsWindow.append(""); if (plotTopN > 0) { // Check the xValues are unique. Since the filters have been sorted by their // numeric value we only need to compare adjacent entries. boolean unique = true; for (int ii = 0; ii < xValues.length - 1; ii++) { if (xValues[ii] == xValues[ii + 1]) { unique = false; break; } } String xAxisName = filterSet.getValueName(); // Check the values all refer to the same property for (Filter filter : filterSet.getFilters()) { if (!xAxisName.equals(filter.getNumericalValueName())) { unique = false; break; } } if (!unique) { // If not unique then renumber them and use an arbitrary label xAxisName = "Filter"; for (int ii = 0; ii < xValues.length; ii++) xValues[ii] = ii + 1; } String title = filterSet.getName(); // Check if a previous filter set had the same name, update if necessary NamedPlot p = getNamedPlot(title); if (p == null) plots.add(new NamedPlot(title, xAxisName, xValues, yValues)); else p.updateValues(xAxisName, xValues, yValues); if (plots.size() > plotTopN) { Collections.sort(plots); p = plots.remove(plots.size() - 1); } } } return count; } private double getCriteria(FractionClassificationResult s) { return getScore(s, criteriaIndex, invertCriteria); } private double getScore(FractionClassificationResult s) { return getScore(s, scoreIndex, invertScore); } private double getScore(FractionClassificationResult s, final int index, final boolean invert) { final double score = getScore(s, index); return (invert) ? -score : score; } private double getScore(FractionClassificationResult s, final int index) { // This order must match the COLUMNS order switch (index) { case 0: return s.getTP() + s.getFP(); case 1: return s.getTP(); case 2: return s.getFP(); case 3: return s.getTN(); case 4: return s.getFN(); case 5: return s.getTPR(); case 6: return s.getTNR(); case 7: return s.getPPV(); case 8: return s.getNPV(); case 9: return s.getFPR(); case 10: return s.getFNR(); case 11: return s.getFDR(); case 12: return s.getAccuracy(); case 13: return s.getMCC(); case 14: return s.getInformedness(); case 15: return s.getMarkedness(); case 16: return s.getRecall(); case 17: return s.getPrecision(); case 18: return s.getFScore(1); case 19: return s.getJaccard(); } return 0; } private boolean requiresInversion(final int index) { switch (index) { case 2: // FP case 4: // FN case 9: // FPR case 10: // FNR case 11: // FDR return true; default: return false; } } private NamedPlot getNamedPlot(String title) { for (NamedPlot p : plots) if (p.name.equals(title)) return p; return null; } private double getMaximum(double[] values) { double max = values[0]; for (int i = 1; i < values.length; i++) { if (values[i] > max) { max = values[i]; } } return max; } private double[] run(Filter filter, List<MemoryPeakResults> resultsList) { FractionClassificationResult r = scoreFilter(filter, resultsList); final double score = getScore(r); final double criteria = getCriteria(r); // Show the result if it achieves the criteria limit if (showResultsTable && criteria >= minCriteria) { String text = createResult(filter, r); if (isHeadless) { IJ.log(text); } else { resultsWindow.append(text); } } return new double[] { score, criteria }; } /** * Score the filter using the results list and the configured fail count * * @param filter * @param resultsList * @return The score */ private FractionClassificationResult scoreFilter(Filter filter, List<MemoryPeakResults> resultsList) { if (failCountRange == 0) return filter.fractionScore(resultsList, failCount); double tp = 0, fp = 0, tn = 0, fn = 0; for (int i = 0; i <= failCountRange; i++) { final FractionClassificationResult r = filter.fractionScore(resultsList, failCount + i); tp += r.getTP(); fp += r.getFP(); tn += r.getTN(); fn += r.getFN(); } return new FractionClassificationResult(tp, fp, tn, fn); } public String createResult(Filter filter, FractionClassificationResult r) { StringBuilder sb = new StringBuilder(resultsPrefix); sb.append(filter.getName()).append(resultsPrefix2).append(resultsPrefix3); int i = 0; add(sb, r.getTP() + r.getFP(), i++); add(sb, r.getTP(), i++); add(sb, r.getFP(), i++); add(sb, r.getTN(), i++); add(sb, r.getFN(), i++); add(sb, r.getTPR(), i++); add(sb, r.getTNR(), i++); add(sb, r.getPPV(), i++); add(sb, r.getNPV(), i++); add(sb, r.getFPR(), i++); add(sb, r.getFNR(), i++); add(sb, r.getFDR(), i++); add(sb, r.getAccuracy(), i++); add(sb, r.getMCC(), i++); add(sb, r.getInformedness(), i++); add(sb, r.getMarkedness(), i++); add(sb, r.getRecall(), i++); add(sb, r.getPrecision(), i++); add(sb, r.getFScore(1), i++); add(sb, r.getJaccard(), i++); // Score relative to the original simulated number of spots // Score the fitting results: // TP are all fit results that can be matched to a spot // FP are all fit results that cannot be matched to a spot // FN = The number of missed spots final double tp = r.getTP(); final double fp = r.getFP(); FractionClassificationResult m = new FractionClassificationResult(tp, fp, 0, simulationParameters.molecules - tp); add(sb, m.getRecall()); add(sb, m.getPrecision()); add(sb, m.getFScore(1)); add(sb, m.getJaccard()); return sb.toString(); } private static void add(StringBuilder sb, String value) { sb.append("\t").append(value); } @SuppressWarnings("unused") private static void add(StringBuilder sb, int value, int i) { if (showColumns[i]) sb.append("\t").append(value); } private static void add(StringBuilder sb, double value, int i) { if (showColumns[i]) add(sb, Utils.rounded(value)); } private static void add(StringBuilder sb, double value) { add(sb, Utils.rounded(value)); } private void saveFilter(Filter filter) { // Save the filter to file String filename = Utils.getFilename("Best_Filter_File", filterFilename); if (filename != null) { List<Filter> filters = new LinkedList<Filter>(); filters.add(filter); FilterSet set = new FilterSet(filter.getName(), filters); List<FilterSet> list = new LinkedList<FilterSet>(); list.add(set); OutputStreamWriter out = null; try { filterFilename = filename; // Append .xml if no suffix if (filename.lastIndexOf('.') < 0) filterFilename += ".xml"; FileOutputStream fos = new FileOutputStream(filterFilename); out = new OutputStreamWriter(fos, "UTF-8"); out.write(XmlUtils.prettyPrintXml(XmlUtils.toXML(list))); } catch (Exception e) { IJ.log("Unable to save the filter sets to file: " + e.getMessage()); } finally { if (out != null) { try { out.close(); } catch (IOException e) { // Ignore } } } } } /** * @param filter */ private void depthAnalysis(Filter filter) { // TODO : This analysis ignores the partial match distance. // Use the score for each result to get a weighted histogram. if (!depthRecallAnalysis || simulationParameters.fixedDepth) return; // Build a histogram of the number of spots at different depths final double[] depths = depthStats.getValues(); final double range = simulationParameters.depth / simulationParameters.a / 2; double[] limits = { -range, range }; final int bins = Math.max(10, simulationParameters.molecules / 50); double[][] h = Utils.calcHistogram(depths, limits[0], limits[1], bins); double[][] h2 = Utils.calcHistogram(depthFitStats.getValues(), limits[0], limits[1], bins); // To get the number of TP at each depth will require that the filter is run // manually to get the results that pass. MemoryPeakResults results = filter.filter(resultsList.get(0), failCount); double[] depths2 = new double[results.size()]; int count = 0; for (PeakResult r : results.getResults()) { if (r.origValue != 0) depths2[count++] = ((DepthPeakResult) r).depth; } depths2 = Arrays.copyOf(depths2, count); // Build a histogram using the same limits double[][] h3 = Utils.calcHistogram(depths2, limits[0], limits[1], bins); // Convert pixel depth to nm for (int i = 0; i < h[0].length; i++) h[0][i] *= simulationParameters.a; limits[0] *= simulationParameters.a; limits[1] *= simulationParameters.a; // Produce a histogram of the number of spots at each depth String title = TITLE + " Depth Histogram"; Plot2 plot = new Plot2(title, "Depth (nm)", "Frequency"); plot.setLimits(limits[0], limits[1], 0, Maths.max(h[1])); plot.setColor(Color.black); plot.addPoints(h[0], h[1], Plot2.BAR); plot.setColor(Color.blue); plot.addPoints(h[0], h2[1], Plot2.BAR); plot.setColor(Color.red); plot.addPoints(h[0], h3[1], Plot2.BAR); plot.addLabel(0, 0, "Black = Spots; Blue = Fitted; Red = Filtered"); PlotWindow pw = Utils.display(title, plot); // Interpolate final double halfBinWidth = (h[0][1] - h[0][0]) * 0.5; // Remove final value of the histogram as this is at the upper limit of the range (i.e. count zero) h[0] = Arrays.copyOf(h[0], h[0].length - 1); h[1] = Arrays.copyOf(h[1], h[0].length); h2[1] = Arrays.copyOf(h2[1], h[0].length); h3[1] = Arrays.copyOf(h3[1], h[0].length); // Use minimum of 3 points for smoothing // Ensure we use at least 15% of data double bandwidth = Math.max(3.0 / h[0].length, 0.15); LoessInterpolator l = new LoessInterpolator(bandwidth, 1); PolynomialSplineFunction spline = l.interpolate(h[0], h[1]); PolynomialSplineFunction spline2 = l.interpolate(h[0], h2[1]); PolynomialSplineFunction spline3 = l.interpolate(h[0], h3[1]); // Increase the number of points to show a smooth curve double[] points = new double[bins * 5]; limits = Maths.limits(h[0]); final double interval = (limits[1] - limits[0]) / (points.length - 1); double[] v = new double[points.length]; double[] v2 = new double[points.length]; double[] v3 = new double[points.length]; for (int i = 0; i < points.length - 1; i++) { points[i] = limits[0] + i * interval; v[i] = spline.value(points[i]); v2[i] = spline2.value(points[i]); v3[i] = spline3.value(points[i]); points[i] += halfBinWidth; } // Final point on the limit of the spline range int ii = points.length - 1; v[ii] = spline.value(limits[1]); v2[ii] = spline2.value(limits[1]); v3[ii] = spline3.value(limits[1]); points[ii] = limits[1] + halfBinWidth; // Calculate recall for (int i = 0; i < v.length; i++) { v2[i] = v2[i] / v[i]; v3[i] = v3[i] / v[i]; } title = TITLE + " Depth Histogram (normalised)"; plot = new Plot2(title, "Depth (nm)", "Recall"); plot.setLimits(limits[0] + halfBinWidth, limits[1] + halfBinWidth, 0, Maths.max(v2)); plot.setColor(Color.blue); plot.addPoints(points, v2, Plot2.LINE); plot.setColor(Color.red); plot.addPoints(points, v3, Plot2.LINE); plot.addLabel(0, 0, "Blue = Fitted; Red = Filtered"); PlotWindow pw2 = Utils.display(title, plot); if (Utils.isNewWindow()) { Point p = pw.getLocation(); p.y += pw.getHeight(); pw2.setLocation(p); } } public class FilterScore implements Comparable<FilterScore> { Filter filter; double score; public FilterScore(Filter filter, double score) { update(filter, score); } public void update(Filter filter, double score) { this.filter = filter; this.score = score; } public int compareTo(FilterScore that) { if (this.score > that.score) return -1; if (this.score < that.score) return 1; return 0; } } public class NamedPlot implements Comparable<NamedPlot> { String name, xAxisName; double[] xValues, yValues; double score; public NamedPlot(String name, String xAxisName, double[] xValues, double[] yValues) { this.name = name; updateValues(xAxisName, xValues, yValues); } public void updateValues(String xAxisName, double[] xValues, double[] yValues) { this.xAxisName = xAxisName; this.xValues = xValues; this.yValues = yValues; this.score = getMaximum(yValues); } public int compareTo(NamedPlot o) { if (score > o.score) return -1; if (score < o.score) return 1; return 0; } } private class DepthPeakResult extends PeakResult { double depth; public DepthPeakResult(int peak, int origX, int origY, float score, double error, float noise, float[] params, double depth) { super(peak, origX, origY, score, error, noise, params, null); this.depth = depth; } } }