Probieren Sie gegebene Punkte stochastisch in einem 3D-Raum mit minimalem Abstand zum nächsten Nachbarn und maximaler Dichte aus

Jan 10 2021

Ich habe nPunkte in einem 3D-Raum. Ich möchte eine Teilmenge von Punkten stochastisch abtasten, wobei alle Abstände zum nächsten Nachbarn größer als sind r. Die Größe der Teilmenge mist unbekannt, aber ich möchte, dass die abgetasteten Punkte so dicht wie möglich sind.

Es gibt ähnliche Fragen, aber es geht nur darum, Punkte zu generieren, anstatt aus bestimmten Punkten zu stichproben.
Generieren Sie zufällige Punkte im 3D-Raum mit minimalem Abstand zum nächsten Nachbarn

3D-Zufallspunkte mit minimalem Abstand zwischen ihnen generieren?

Angenommen, ich habe 300 zufällige 3D-Punkte.

import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))

Was ist der schnellste Weg, um eine Teilmenge von mPunkten mit minimalem Abstand zum nächsten Nachbarn zu erhalten und r = 1gleichzeitig zu maximieren m?

Antworten

2 DavidEisenstat Jan 14 2021 at 08:00

Es gibt wahrscheinlich ein effizientes Bicriteria-Approximationsschema, aber warum sollte man sich die Mühe machen, wenn die Ganzzahlprogrammierung im Durchschnitt so schnell ist?

import numpy as np

n = 300
points = np.random.uniform(0, 10, size=(n, 3))

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver("SCIP")
variables = [solver.BoolVar("x[{}]".format(i)) for i in range(n)]
solver.Maximize(sum(variables))
for j, q in enumerate(points):
    for i, p in enumerate(points[:j]):
        if np.linalg.norm(p - q) <= 1:
            solver.Add(variables[i] + variables[j] <= 1)
solver.EnableOutput()
solver.Solve()
print(len([i for (i, variable) in enumerate(variables) if variable.SolutionValue()]))
1 DanielF Jan 18 2021 at 16:16

Dies ist nicht optimal groß für eine Teilmenge, sollte aber nahe sein, ohne sehr lange zu dauern KDTree, um die Entfernungsberechnungen zu optimieren:

from scipy.spatial import KDTree
import numpy as np

def space_sample(n = 300, low = 0, high = 10, dist = 1):
    points = np.random.uniform(low, high, size=(n, 3))
    k = KDTree(points)
    pairs = np.array(list(k.query_pairs(dist)))
    
    def reduce_pairs(pairs, remove = []):  #iteratively remove the most connected node
        p = pairs[~np.isin(pairs, remove).any(1)]
        if p.size >0:
            count = np.bincount(p.flatten(), minlength = n)
            r = remove + [count.argmax()]
            return reduce_pairs(p, r)
        else:
            return remove
    
    return np.array([p for i, p in enumerate(points) if not(i in reduce_pairs(pairs))])

subset = space_sample()

Das iterative Entfernen des am meisten verbundenen Knotens ist nicht optimal (behält ungefähr 200 der 300 Punkte bei), ist jedoch wahrscheinlich der schnellste Algorithmus, der nahezu optimal ist (das einzige, was schneller ist, ist das zufällige Entfernen). Sie könnten @njit reduce_pairsdiesen Teil möglicherweise schneller machen (werde es versuchen, wenn ich später Zeit habe).

ShaunHan Jan 19 2021 at 07:07

Testen der Antwort von @David Eisenstat mit 30 gegebenen Punkten:

from ortools.linear_solver import pywraplp
import numpy as np

def subset_David_Eisenstat(points, r):
    solver = pywraplp.Solver.CreateSolver("SCIP")
    variables = [solver.BoolVar("x[{}]".format(i)) for i in range(len(points))]
    solver.Maximize(sum(variables))
    for j, q in enumerate(points):
        for i, p in enumerate(points[:j]):
            if np.linalg.norm(p - q) <= r:
                solver.Add(variables[i] + variables[j] <= 1)
    solver.EnableOutput()
    solver.Solve()
    indices = [i for (i, variable) in enumerate(variables) if variable.SolutionValue()]
    return points[indices]

points = np.array(
[[ 7.32837882, 12.12765786, 15.01412241],
 [ 8.25164031, 11.14830379, 15.01412241],
 [ 8.21790113, 13.05647987, 13.05647987],
 [ 7.30031002, 13.08276009, 14.05452502],
 [ 9.18056467, 12.0800735 , 13.05183844],
 [ 9.17929647, 11.11270337, 14.03027534],
 [ 7.64737905, 11.48906945, 15.34274827],
 [ 7.01315886, 12.77870699, 14.70301668],
 [ 8.88132414, 10.81243313, 14.68685022],
 [ 7.60617372, 13.39792166, 13.39792166],
 [ 8.85967682, 12.72946394, 12.72946394],
 [ 9.50060705, 11.43361294, 13.37866092],
 [ 8.21790113, 12.08471494, 14.02824481],
 [ 7.32837882, 12.12765786, 16.98587759],
 [ 8.25164031, 11.14830379, 16.98587759],
 [ 7.30031002, 13.08276009, 17.94547498],
 [ 8.21790113, 13.05647987, 18.94352013],
 [ 9.17929647, 11.11270337, 17.96972466],
 [ 9.18056467, 12.0800735 , 18.94816156],
 [ 7.64737905, 11.48906945, 16.65725173],
 [ 7.01315886, 12.77870699, 17.29698332],
 [ 8.88132414, 10.81243313, 17.31314978],
 [ 7.60617372, 13.39792166, 18.60207834],
 [ 8.85967682, 12.72946394, 19.27053606],
 [ 9.50060705, 11.43361294, 18.62133908],
 [ 8.21790113, 12.08471494, 17.97175519],
 [ 7.32837882, 15.01412241, 12.12765786],
 [ 8.25164031, 15.01412241, 11.14830379],
 [ 7.30031002, 14.05452502, 13.08276009],
 [ 9.18056467, 13.05183844, 12.0800735 ],])

Wenn der erwartete Mindestabstand 1 beträgt:

subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18

Überprüfen Sie den Mindestabstand:

from scipy.spatial.distance import cdist
dist = cdist(subset1, subset1, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 1.3285513450926985

Ändern Sie den erwarteten Mindestabstand auf 2:

subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10

Überprüfen Sie den Mindestabstand:

from scipy.spatial.distance import cdist
dist = cdist(subset2, subset2, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 2.0612041004376223