using StatsBase

"""
    smd_combined = combinelocalizations(smd::SMD)

Combine clustered localizations in `smd` into higher precision localizations.
"""
function combinelocalizations(smd::SMD)
    # Isolate and sort some arrays from `smd`.
    connectID = smd.connectID
    sortindices = sortperm(connectID)
    connectID = connectID[sortindices]
    x = smd.x[sortindices]
    y = smd.y[sortindices]
    x_se = smd.x_se[sortindices]
    y_se = smd.y_se[sortindices]
    framenum = smd.framenum[sortindices]
    datasetnum = smd.datasetnum[sortindices]

    # Loop over clusters and combine their localization data as appropriate.
    nperID = counts(connectID)
    nperID = nperID[nperID .!= 0]
    ncumulative = [0; cumsum(nperID)]
    nclusters = length(nperID)
    smd_combined = deepcopy(smd)
    for nn = 1:nclusters
        indices = (1:nperID[nn]) .+ ncumulative[nn]
        smd_combined.x[nn] = StatsBase.mean(x[indices], weights(1.0/x_se[indices].^2))
        smd_combined.y[nn] = StatsBase.mean(y[indices], weights(1.0/y_se[indices].^2))
        smd_combined.x_se[nn] = sqrt(1.0 / sum(1.0./x_se[indices].^2))
        smd_combined.y_se[nn] = sqrt(1.0 / sum(1.0./y_se[indices].^2))
        smd_combined.connectID[nn] = connectID[indices[1]]
        smd_combined.framenum[nn] = framenum[indices[1]]
        smd_combined.datasetnum[nn] = datasetnum[indices[1]]
    end
    smd_combined = isolatesmd(smd_combined, 1:nclusters)

    return smd_combined
end