gdsc.core.match.MatchCalculator.java Source code

Java tutorial

Introduction

Here is the source code for gdsc.core.match.MatchCalculator.java

Source

package gdsc.core.match;

/*----------------------------------------------------------------------------- 
 * 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 java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
import java.util.List;

import org.apache.commons.math3.util.FastMath;

/**
 * Calculates the match between a set of predicted points and the actual points.
 */
public class MatchCalculator {
    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @return The match results
     */
    public static MatchResult analyseResults2D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold) {
        return analyseResults2D(actualPoints, predictedPoints, dThreshold, null, null, null, null);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @return The match results
     */
    public static MatchResult analyseResults2D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN) {
        return analyseResults2D(actualPoints, predictedPoints, dThreshold, TP, FP, FN, null);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @param matches
     *            The matched true positives (point1 = actual, point2 = predicted)
     * @return The match results
     */
    public static MatchResult analyseResults2D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN,
            List<PointPair> matches) {
        final int predictedPointsLength = (predictedPoints != null) ? predictedPoints.length : 0;
        final int actualPointsLength = (actualPoints != null) ? actualPoints.length : 0;

        // If the number of possible pairs is small then use a one pass algorithm
        if (predictedPointsLength * actualPointsLength < 100000)
            return analyseResults2DSinglePass(actualPoints, predictedPoints, dThreshold, TP, FP, FN, matches);

        dThreshold *= dThreshold; // We will use the squared distance

        int tp = 0; // true positives (actual with matched predicted point)
        int fp = predictedPointsLength; // false positives (actual with no matched predicted point)
        int fn = actualPointsLength; // false negatives (predicted point with no actual point)
        double rmsd = 0;

        clear(TP, FP, FN, matches);
        if (predictedPointsLength == 0 || actualPointsLength == 0) {
            if (FP != null)
                FP.addAll(asList(predictedPoints));
            if (FN != null)
                FN.addAll(asList(actualPoints));
            return new MatchResult(tp, fp, fn, rmsd);
        }

        // loop over the two arrays assigning the closest unassigned pair
        final boolean[] resultAssignment = new boolean[predictedPointsLength];
        final boolean[] roiAssignment = new boolean[fn];
        final ArrayList<ImmutableAssignment> assignments = new ArrayList<ImmutableAssignment>(
                predictedPointsLength);

        int[] falsePositives = null, falseNegatives = null;
        if (FP != null) {
            falsePositives = ascendingArray(predictedPointsLength);
        }
        if (FN != null) {
            falseNegatives = ascendingArray(actualPointsLength);
        }

        do {
            assignments.clear();

            // Process each result
            for (int predictedId = predictedPointsLength; predictedId-- > 0;) {
                if (resultAssignment[predictedId])
                    continue; // Already assigned

                final float x = predictedPoints[predictedId].getX();
                final float y = predictedPoints[predictedId].getY();

                // Find closest ROI point
                float d2Min = (float) dThreshold; //Float.MAX_VALUE;
                int targetId = -1;
                for (int actualId = actualPointsLength; actualId-- > 0;) {
                    if (roiAssignment[actualId])
                        continue; // Already assigned

                    Coordinate actualPoint = actualPoints[actualId];

                    // Calculate in steps for increased speed (allows early exit)
                    float dx = actualPoint.getX() - x;
                    dx *= dx;
                    if (dx <= d2Min) {
                        float dy = actualPoint.getY() - y;
                        dy *= dy;
                        if (dy <= d2Min) {
                            final float d2 = dx + dy;
                            if (d2 <= d2Min) {
                                d2Min = d2;
                                targetId = actualId;
                            }
                        }
                    }
                }

                // Store closest ROI point
                if (targetId > -1) {
                    assignments.add(new ImmutableAssignment(targetId, predictedId, d2Min));
                }
            }

            // If there are assignments
            if (!assignments.isEmpty()) {
                // Pick the closest pair to be assigned
                AssignmentComparator.sort(assignments);

                // Process in order
                for (ImmutableAssignment closest : assignments) {
                    // Skip those that have already been assigned since this will be a lower score.
                    // Note at least one assignment should be processed as potential assignments are made 
                    // using only unassigned points.
                    if (resultAssignment[closest.getPredictedId()] || roiAssignment[closest.getTargetId()])
                        continue;
                    resultAssignment[closest.getPredictedId()] = true;
                    roiAssignment[closest.getTargetId()] = true;

                    // If within accuracy then classify as a match
                    if (closest.getDistance() <= dThreshold) {
                        tp++;
                        fn--;
                        fp--;
                        rmsd += closest.getDistance(); // Already a squared distance

                        if (TP != null) {
                            TP.add(predictedPoints[closest.getPredictedId()]);
                        }
                        if (FP != null) {
                            falsePositives[closest.getPredictedId()] = -1;
                        }
                        if (FN != null) {
                            falseNegatives[closest.getTargetId()] = -1;
                        }
                        if (matches != null)
                            matches.add(new PointPair(actualPoints[closest.getTargetId()],
                                    predictedPoints[closest.getPredictedId()]));
                    } else {
                        // No more assignments within the distance threshold
                        break;
                    }
                }
            }

        } while (!assignments.isEmpty());

        // Add to lists
        if (FP != null) {
            for (int i = 0; i < predictedPointsLength; i++) {
                if (falsePositives[i] >= 0)
                    FP.add(predictedPoints[i]);
            }
        }
        if (FN != null) {
            for (int i = 0; i < actualPointsLength; i++) {
                if (falseNegatives[i] >= 0)
                    FN.add(actualPoints[i]);
            }
        }

        if (tp > 0)
            rmsd = Math.sqrt(rmsd / tp);
        return new MatchResult(tp, fp, fn, rmsd);
    }

    private static <T, U, V> void clear(List<T> TP, List<T> FP, List<U> FN, List<V> matches) {
        clear(TP);
        clear(FP);
        clear(FN);
        clear(matches);
    }

    private static <T> void clear(List<T> list) {
        if (list != null)
            list.clear();
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * <p>
     * Use a single pass algorithm suitable if the total number of possible pairs is small (<100000)
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @param matches
     *            The matched true positives (point1 = actual, point2 = predicted)
     * @return The match results
     */
    public static MatchResult analyseResults2DSinglePass(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN,
            List<PointPair> matches) {
        dThreshold *= dThreshold; // We will use the squared distance

        final int predictedPointsLength = (predictedPoints != null) ? predictedPoints.length : 0;
        final int actualPointsLength = (actualPoints != null) ? actualPoints.length : 0;

        int tp = 0; // true positives (actual with matched predicted point)
        int fp = predictedPointsLength; // false positives (actual with no matched predicted point)
        int fn = actualPointsLength; // false negatives (predicted point with no actual point)
        double rmsd = 0;

        clear(TP, FP, FN, matches);
        if (predictedPointsLength == 0 || actualPointsLength == 0) {
            if (FP != null)
                FP.addAll(asList(predictedPoints));
            if (FN != null)
                FN.addAll(asList(actualPoints));
            return new MatchResult(tp, fp, fn, rmsd);
        }

        // loop over the two arrays assigning the closest unassigned pair
        final ArrayList<ImmutableAssignment> assignments = new ArrayList<ImmutableAssignment>(
                predictedPointsLength);

        for (int predictedId = predictedPointsLength; predictedId-- > 0;) {
            final float x = predictedPoints[predictedId].getX();
            final float y = predictedPoints[predictedId].getY();
            for (int actualId = actualPointsLength; actualId-- > 0;) {
                final double d2 = actualPoints[actualId].distance2(x, y);
                if (d2 <= dThreshold) {
                    assignments.add(new ImmutableAssignment(actualId, predictedId, d2));
                }
            }
        }

        AssignmentComparator.sort(assignments);

        final boolean[] predictedAssignment = new boolean[predictedPointsLength];
        final boolean[] actualAssignment = new boolean[actualPointsLength];

        for (ImmutableAssignment a : assignments) {
            if (!actualAssignment[a.getTargetId()]) {
                if (!predictedAssignment[a.getPredictedId()]) {
                    actualAssignment[a.getTargetId()] = true;
                    predictedAssignment[a.getPredictedId()] = true;
                    tp++;
                    fn--;
                    fp--;
                    rmsd += a.getDistance(); // Already a squared distance
                    if (matches != null)
                        matches.add(
                                new PointPair(actualPoints[a.getTargetId()], predictedPoints[a.getPredictedId()]));
                    if (TP != null)
                        TP.add(predictedPoints[a.getPredictedId()]);
                }
            }
        }

        // Add to lists
        if (FP != null) {
            for (int i = 0; i < predictedPointsLength; i++) {
                if (!predictedAssignment[i])
                    FP.add(predictedPoints[i]);
            }
        }
        if (FN != null) {
            for (int i = 0; i < actualPointsLength; i++) {
                if (!actualAssignment[i])
                    FN.add(actualPoints[i]);
            }
        }

        if (tp > 0)
            rmsd = Math.sqrt(rmsd / tp);
        return new MatchResult(tp, fp, fn, rmsd);
    }

    private static Collection<Coordinate> asList(Coordinate[] points) {
        if (points != null)
            return Arrays.asList(points);
        return new ArrayList<Coordinate>(0);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @return The match results
     */
    public static MatchResult analyseResults3D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold) {
        return analyseResults3D(actualPoints, predictedPoints, dThreshold, null, null, null, null);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @return The match results
     */
    public static MatchResult analyseResults3D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN) {
        return analyseResults3D(actualPoints, predictedPoints, dThreshold, TP, FP, FN, null);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @param matches
     *            The matched true positives (point1 = actual, point2 = predicted)
     * @return The match results
     */
    public static MatchResult analyseResults3D(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN,
            List<PointPair> matches) {
        final int predictedPointsLength = (predictedPoints != null) ? predictedPoints.length : 0;
        final int actualPointsLength = (actualPoints != null) ? actualPoints.length : 0;

        // If the number of possible pairs is small then use a one pass algorithm
        if (predictedPointsLength * actualPointsLength < 100000)
            return analyseResults3DSinglePass(actualPoints, predictedPoints, dThreshold, TP, FP, FN, matches);

        dThreshold *= dThreshold; // We will use the squared distance

        int tp = 0; // true positives (actual with matched predicted point)
        int fp = predictedPointsLength; // false positives (actual with no matched predicted point)
        int fn = actualPointsLength; // false negatives (predicted point with no actual point)
        double rmsd = 0;

        clear(TP, FP, FN, matches);
        if (predictedPointsLength == 0 || actualPointsLength == 0) {
            if (FP != null)
                FP.addAll(asList(predictedPoints));
            if (FN != null)
                FN.addAll(asList(actualPoints));
            return new MatchResult(tp, fp, fn, rmsd);
        }

        // loop over the two arrays assigning the closest unassigned pair
        final boolean[] resultAssignment = new boolean[predictedPointsLength];
        final boolean[] roiAssignment = new boolean[fn];
        final ArrayList<ImmutableAssignment> assignments = new ArrayList<ImmutableAssignment>(
                predictedPointsLength);

        int[] falsePositives = null, falseNegatives = null;
        if (FP != null) {
            falsePositives = ascendingArray(predictedPointsLength);
        }
        if (FN != null) {
            falseNegatives = ascendingArray(actualPointsLength);
        }

        do {
            assignments.clear();

            // Process each result
            for (int predictedId = predictedPointsLength; predictedId-- > 0;) {
                if (resultAssignment[predictedId])
                    continue; // Already assigned

                final float x = predictedPoints[predictedId].getX();
                final float y = predictedPoints[predictedId].getY();
                final float z = predictedPoints[predictedId].getZ();

                // Find closest ROI point
                float d2Min = (float) dThreshold; //Float.MAX_VALUE;
                int targetId = -1;
                for (int actualId = actualPointsLength; actualId-- > 0;) {
                    if (roiAssignment[actualId])
                        continue; // Already assigned

                    Coordinate actualPoint = actualPoints[actualId];

                    // Calculate in steps for increased speed (allows early exit)
                    float dx = actualPoint.getX() - x;
                    dx *= dx;
                    if (dx <= d2Min) {
                        float dy = actualPoint.getY() - y;
                        dy *= dy;
                        if (dy <= d2Min) {
                            float dz = actualPoint.getZ() - z;
                            dz *= dz;
                            if (dz <= d2Min) {
                                final float d2 = dx + dy + dz;
                                if (d2 <= d2Min) {
                                    d2Min = d2;
                                    targetId = actualId;
                                }
                            }
                        }
                    }
                }

                // Store closest ROI point
                if (targetId > -1) {
                    assignments.add(new ImmutableAssignment(targetId, predictedId, d2Min));
                }
            }

            // If there are assignments
            if (!assignments.isEmpty()) {
                // Pick the closest pair to be assigned
                AssignmentComparator.sort(assignments);

                // Process in order
                for (ImmutableAssignment closest : assignments) {
                    // Skip those that have already been assigned since this will be a lower score.
                    // Note at least one assignment should be processed as potential assignments are made 
                    // using only unassigned points.
                    if (resultAssignment[closest.getPredictedId()] || roiAssignment[closest.getTargetId()])
                        continue;

                    resultAssignment[closest.getPredictedId()] = true;
                    roiAssignment[closest.getTargetId()] = true;

                    // If within accuracy then classify as a match
                    if (closest.getDistance() <= dThreshold) {
                        tp++;
                        fn--;
                        fp--;
                        rmsd += closest.getDistance();

                        if (TP != null) {
                            TP.add(predictedPoints[closest.getPredictedId()]);
                        }
                        if (FP != null) {
                            falsePositives[closest.getPredictedId()] = -1;
                        }
                        if (FN != null) {
                            falseNegatives[closest.getTargetId()] = -1;
                        }
                        if (matches != null)
                            matches.add(new PointPair(actualPoints[closest.getTargetId()],
                                    predictedPoints[closest.getPredictedId()]));
                    } else {
                        // No more assignments within the distance threshold
                        break;
                    }
                }
            }

        } while (!assignments.isEmpty());

        // Add to lists
        if (FP != null) {
            for (int i = 0; i < predictedPointsLength; i++) {
                if (falsePositives[i] >= 0)
                    FP.add(predictedPoints[i]);
            }
        }
        if (FN != null) {
            for (int i = 0; i < actualPointsLength; i++) {
                if (falseNegatives[i] >= 0)
                    FN.add(actualPoints[i]);
            }
        }

        if (tp > 0)
            rmsd = Math.sqrt(rmsd / tp);
        return new MatchResult(tp, fp, fn, rmsd);
    }

    /**
     * Calculate the match results for the given actual and predicted points.
     * Points that are within the distance threshold are identified as a match.
     * The number of true positives, false positives and false negatives are calculated.
     * <p>
     * Use a single pass algorithm suitable if the total number of possible pairs is small (<100000)
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @param matches
     *            The matched true positives (point1 = actual, point2 = predicted)
     * @return The match results
     */
    public static MatchResult analyseResults3DSinglePass(Coordinate[] actualPoints, Coordinate[] predictedPoints,
            double dThreshold, List<Coordinate> TP, List<Coordinate> FP, List<Coordinate> FN,
            List<PointPair> matches) {
        dThreshold *= dThreshold; // We will use the squared distance

        final int predictedPointsLength = (predictedPoints != null) ? predictedPoints.length : 0;
        final int actualPointsLength = (actualPoints != null) ? actualPoints.length : 0;

        int tp = 0; // true positives (actual with matched predicted point)
        int fp = predictedPointsLength; // false positives (actual with no matched predicted point)
        int fn = actualPointsLength; // false negatives (predicted point with no actual point)
        double rmsd = 0;

        clear(TP, FP, FN, matches);
        if (predictedPointsLength == 0 || actualPointsLength == 0) {
            if (FP != null)
                FP.addAll(asList(predictedPoints));
            if (FN != null)
                FN.addAll(asList(actualPoints));
            return new MatchResult(tp, fp, fn, rmsd);
        }

        // loop over the two arrays assigning the closest unassigned pair
        final ArrayList<ImmutableAssignment> assignments = new ArrayList<ImmutableAssignment>(
                predictedPointsLength);

        for (int predictedId = predictedPointsLength; predictedId-- > 0;) {
            final float x = predictedPoints[predictedId].getX();
            final float y = predictedPoints[predictedId].getY();
            final float z = predictedPoints[predictedId].getZ();
            for (int actualId = actualPointsLength; actualId-- > 0;) {
                final double d2 = actualPoints[actualId].distance2(x, y, z);
                if (d2 <= dThreshold) {
                    assignments.add(new ImmutableAssignment(actualId, predictedId, d2));
                }
            }
        }

        AssignmentComparator.sort(assignments);

        final boolean[] predictedAssignment = new boolean[predictedPointsLength];
        final boolean[] actualAssignment = new boolean[actualPointsLength];

        for (ImmutableAssignment a : assignments) {
            if (!actualAssignment[a.getTargetId()]) {
                if (!predictedAssignment[a.getPredictedId()]) {
                    actualAssignment[a.getTargetId()] = true;
                    predictedAssignment[a.getPredictedId()] = true;
                    tp++;
                    fn--;
                    fp--;
                    rmsd += a.getDistance(); // Already a squared distance
                    if (matches != null)
                        matches.add(
                                new PointPair(actualPoints[a.getTargetId()], predictedPoints[a.getPredictedId()]));
                    if (TP != null)
                        TP.add(predictedPoints[a.getPredictedId()]);
                }
            }
        }

        // Add to lists
        if (FP != null) {
            for (int i = 0; i < predictedPointsLength; i++) {
                if (!predictedAssignment[i])
                    FP.add(predictedPoints[i]);
            }
        }
        if (FN != null) {
            for (int i = 0; i < actualPointsLength; i++) {
                if (!actualAssignment[i])
                    FN.add(actualPoints[i]);
            }
        }

        if (tp > 0)
            rmsd = Math.sqrt(rmsd / tp);
        return new MatchResult(tp, fp, fn, rmsd);
    }

    private static int[] ascendingArray(int length) {
        int[] array = new int[length];
        for (int i = 0; i < length; i++)
            array[i] = i;
        return array;
    }

    /**
     * Calculate the match results for the given actual and predicted fluorophore pulses.
     * Points that are within the distance threshold are identified as a match. The score is calculated using half the
     * distance threshold and the overlap in time. Assignments are made using the highest scoring matches.
     * <p>
     * The total score is stored in the RMSD field of the MatchResult. The number of true positives, false positives and
     * false negatives are calculated.
     * 
     * @param actualPoints
     * @param predictedPoints
     * @param dThreshold
     *            The distance threshold
     * @param TP
     *            True Positives
     * @param FP
     *            False Positives
     * @param FN
     *            False Negatives
     * @param matches
     *            The matched true positives (point1 = actual, point2 = predicted)
     * @return The match results
     */
    public static MatchResult analyseResults2D(Pulse[] actualPoints, Pulse[] predictedPoints,
            final double dThreshold, List<Pulse> TP, List<Pulse> FP, List<Pulse> FN, List<PointPair> matches) {
        // We will use the squared distance for speed
        final double halfDThreshold2 = (dThreshold * 0.5) * (dThreshold * 0.5);
        final float floatDThreshold = (float) dThreshold;
        final double dThreshold2 = dThreshold * dThreshold;

        final int predictedPointsLength = (predictedPoints != null) ? predictedPoints.length : 0;
        final int actualPointsLength = (actualPoints != null) ? actualPoints.length : 0;

        int tp = 0; // true positives (actual with matched predicted point)
        int fp = predictedPointsLength; // false positives (actual with no matched predicted point)
        int fn = actualPointsLength; // false negatives (predicted point with no actual point)
        double score = 0;

        clear(TP, FP, FN, matches);
        if (predictedPointsLength == 0 || actualPointsLength == 0) {
            if (FP != null)
                FP.addAll(asList(predictedPoints));
            if (FN != null)
                FN.addAll(asList(actualPoints));
            return new MatchResult(tp, fp, fn, score);
        }

        // loop over the two arrays assigning the closest unassigned pair
        boolean[] resultAssignment = new boolean[predictedPointsLength];
        boolean[] roiAssignment = new boolean[fn];
        ArrayList<ImmutableAssignment> assignments = new ArrayList<ImmutableAssignment>(predictedPointsLength);

        int[] falsePositives = null, falseNegatives = null;
        if (FP != null) {
            falsePositives = ascendingArray(predictedPointsLength);
        }
        if (FN != null) {
            falseNegatives = ascendingArray(actualPointsLength);
        }

        // Sort by time to allow efficient looping
        Arrays.sort(actualPoints);
        Arrays.sort(predictedPoints);

        // Pre-calculate all-vs-all distance matrix if it can fit in memory
        int size = predictedPointsLength * actualPointsLength;
        float[][] dMatrix = null;
        if (size < 200 * 200) {
            dMatrix = new float[predictedPointsLength][actualPointsLength];
            for (int predictedId = 0; predictedId < predictedPointsLength; predictedId++) {
                final float x = predictedPoints[predictedId].getX();
                final float y = predictedPoints[predictedId].getY();
                for (int actualId = 0; actualId < actualPointsLength; actualId++) {
                    dMatrix[predictedId][actualId] = (float) actualPoints[actualId].distance2(x, y);
                }
            }
        }

        do {
            assignments.clear();

            // Process each result
            for (int predictedId = 0; predictedId < predictedPointsLength; predictedId++) {
                if (resultAssignment[predictedId])
                    continue; // Already assigned

                final float x = predictedPoints[predictedId].getX();
                final float y = predictedPoints[predictedId].getY();
                final int start = predictedPoints[predictedId].getStart();
                final int end = predictedPoints[predictedId].getEnd();

                // Find first overlapping pulse
                int actualId = 0;
                while (actualId < actualPointsLength && actualPoints[actualId].getEnd() < start) {
                    actualId++;
                }

                // Find highest scoring point within the distance limit
                double scoreMax = 0;
                int targetId = -1;
                for (; actualId < actualPointsLength; actualId++) {
                    if (roiAssignment[actualId])
                        continue; // Already assigned
                    if (actualPoints[actualId].getStart() > end)
                        break; // No more overlap in time

                    double d2;
                    if (dMatrix != null) {
                        d2 = dMatrix[predictedId][actualId];
                    } else {
                        Coordinate actualPoint = actualPoints[actualId];

                        // Calculate in steps for increased speed (allows early exit)
                        float dx = abs(actualPoint.getX() - x);
                        if (dx > floatDThreshold)
                            continue;
                        float dy = abs(actualPoint.getY() - y);
                        if (dy > floatDThreshold)
                            continue;
                        d2 = dx * dx + dy * dy;
                    }

                    // Do we need to exclude using the distance threshold? This is useful for binary classification
                    // but will truncate the continuous nature of the score.
                    if (d2 > dThreshold2)
                        continue;

                    double s = predictedPoints[predictedId].score(actualPoints[actualId], d2, halfDThreshold2);
                    if (scoreMax < s) {
                        scoreMax = s;
                        targetId = actualId;
                    }
                }

                // Store highest scoring point
                if (targetId > -1) {
                    assignments.add(new ImmutableAssignment(targetId, predictedId, scoreMax));
                }
            }

            // If there are assignments
            if (!assignments.isEmpty()) {
                // Process highest scoring first
                AssignmentComparator.sort(assignments);
                Collections.reverse(assignments);

                // Process in order of score
                for (ImmutableAssignment closest : assignments) {
                    // Skip those that have already been assigned since this will be a lower score.
                    // Note at least one assignment should be processed as potential assignments are made 
                    // using only unassigned points.
                    if (resultAssignment[closest.getPredictedId()] || roiAssignment[closest.getTargetId()])
                        continue;

                    resultAssignment[closest.getPredictedId()] = true;
                    roiAssignment[closest.getTargetId()] = true;

                    tp++;
                    fn--;
                    fp--;
                    score += closest.getDistance(); // This is the scoreMax (not the distance)

                    if (TP != null) {
                        TP.add(predictedPoints[closest.getPredictedId()]);
                    }
                    if (FP != null) {
                        falsePositives[closest.getPredictedId()] = -1;
                    }
                    if (FN != null) {
                        falseNegatives[closest.getTargetId()] = -1;
                    }
                    if (matches != null)
                        matches.add(new PointPair(actualPoints[closest.getTargetId()],
                                predictedPoints[closest.getPredictedId()]));
                }
            }

        } while (!assignments.isEmpty());

        // Add to lists
        if (FP != null) {
            for (int i = 0; i < predictedPointsLength; i++) {
                if (falsePositives[i] >= 0)
                    FP.add(predictedPoints[i]);
            }
        }
        if (FN != null) {
            for (int i = 0; i < actualPointsLength; i++) {
                if (falseNegatives[i] >= 0)
                    FN.add(actualPoints[i]);
            }
        }

        // Every time-point has the chance to contribute to the score.
        // Normalise score by the maximum of the number of actual/predicted time points.
        // This penalises too few or too many predictions
        int p1 = countTimePoints(actualPoints);
        int p2 = countTimePoints(predictedPoints);
        score /= FastMath.max(p1, p2);

        return new MatchResult(tp, fp, fn, score);
    }

    private static int countTimePoints(Pulse[] actualPoints) {
        int p1 = 0;
        for (Pulse p : actualPoints)
            p1 += p.getEnd() - p.getStart() + 1;
        return p1;
    }

    private static Collection<Pulse> asList(Pulse[] points) {
        if (points != null)
            return Arrays.asList(points);
        return new ArrayList<Pulse>(0);
    }

    private static float abs(final float f) {
        return (f < 0) ? -f : f;
    }
}