#алгоритм #любой_язык #кластеризация
Описание задачи
Определить геометрический центр группы 3D-объектов для её использования в качестве
опорной точки (pivot point). При этом в модели присутствует основная группа объектов,
расположеная в непосредственной близости друг от друга, и "лишние" объекты, сильно
удаленные от этой группы.
Если использовать простой способ определения геометрического центра всех объектов,
то получаемая координата может быть сильно удалена от основной группы объектов (псевдокод):
double cx = 0, cy = 0, cz = 0;
for (auto obj : objects) {
cx += obj->x;
cy += obj->y;
cz += obj->z;
}
cx /= objects.size();
cy /= objects.size();
cz /= objects.size();
С помощью какого алгоритма можно найти основную группу близко расположенных объектов,
отфильтровав все остальные далеко расположенные одиночные объекты, либо небольшие группы
объектов?
Тестовый набор данных
Далее привожу простой пример решаемой задачи.
Визуальное представление (синим цветом выделены объекты, которые необходимо отфильтровать):
Координаты объектов:
--------------------------------------------------------------------
| X | Z | Y |
--------------------------------------------------------------------
| 2.709080934524536 | -1.9362084865570068 | -1.0364959239959717 |
--------------------------------------------------------------------
| 2.923645257949829 | 1.1126995086669922 | 1.5516517162322998 |
--------------------------------------------------------------------
| 0.5295989513397217 | 3.0490918159484863 | 0.1415269374847412 |
--------------------------------------------------------------------
| -1.4740893840789795 | 1.233896255493164 | -3.516988515853882 |
--------------------------------------------------------------------
| -2.124799966812134 | 0.3458544611930847 | -4.8596062660217285 |
--------------------------------------------------------------------
| -1.0827465057373047 | -0.8037613034248352 | -4.489010810852051 |
--------------------------------------------------------------------
| -0.7635605335235596 | 0.03742170333862305 | -3.5328845977783203 |
--------------------------------------------------------------------
| -1.1618667840957642 | -0.4031369686126709 | -3.6103110313415527 |
--------------------------------------------------------------------
| -1.5830864906311035 | 0.296390175819397 | -3.8527112007141113 |
--------------------------------------------------------------------
| -1.7283073663711548 | -0.2046224325895309 | -4.165205001831055 |
--------------------------------------------------------------------
Ответы
Ответ 1
Вы пытаетесь получить среднее арифметическое значение для центра тяжести вместо средневзвешенного: c = sum(ciVi) / sum Vi. т.е. не учитываете вес (или объём) объектов objects->weight() в первой сумме. Правильный вариант: double cx = 0, cy = 0, cz = 0, cw = 0; for (auto obj : objects) { cx += obj->x*obj->weight; cy += obj->y*obj->weight; cz += obj->z*obj->weight; cw += obj->weight; } cx /= cw; cy /= cw; cz /= cw; Если свойства objects->weight у объектов нет, его следует ввести. Если ввести такое свойство невозможно, то робастной (грубой, но устойчивой) оценкой является медиана. Т.е. надо отсортировать массив (независимо по каждой из координат) и взять средние по ранжиру (если количество элементов нечётное) или полусуммы двух средних (если количество элементов чётное). А в рамках трёхмерной схемы нужно выбирать такой объект, сумма расстояний до которого от остальных объектов минимальна. При этом возникает выбор, что суммировать: расстояния или их квадрат? Это определяется параметрами разброса объектов. если его распределение близко к нормальному, то лучше суммировать квадраты. Если же оно скорее хаотическое (импульсное) или статистика мала, то лучше суммировать расстояния (робастная оценка). Программа с суммированием расстояний и независимыми медианами: $xzy = [ [ 2.709080934524536 , -1.9362084865570068 , -1.0364959239959717 ], [ 2.923645257949829 , 1.1126995086669922 , 1.5516517162322998 ], [ 0.5295989513397217 , 3.0490918159484863 , 0.1415269374847412 ], [ -1.4740893840789795 , 1.233896255493164 , -3.516988515853882 ], [ -2.124799966812134 , 0.3458544611930847 , -4.8596062660217285 ], [ -1.0827465057373047 , -0.8037613034248352 , -4.489010810852051 ], [ -0.7635605335235596 , 0.03742170333862305 , -3.5328845977783203 ], [ -1.1618667840957642 , -0.4031369686126709 , -3.6103110313415527 ], [ -1.5830864906311035 , 0.296390175819397 , -3.8527112007141113 ], [ -1.7283073663711548 , -0.2046224325895309 , -4.165205001831055 ] ]; function min_distances($arr){ $sum_min = 999999.9; $key_min = -1; //var_dump($arr); $brr = $arr; foreach($arr as $key=>$base){ $base_x = $base[0]; $base_y = $base[2]; $base_z = $base[1]; $sum_dist = 0; foreach($arr as $another){ $dx = $another[0]-$base[0]; $dz = $another[1]-$base[1]; $dy = $another[2]-$base[2]; $sum_dist += sqrt($dx*$dx+$dy*$dy+$dz*$dz); } if($sum_dist < $sum_min){ $sum_min = $sum_dist; $key_min = $key; } } return $key_min; } $key_min = min_distances($xzy); printf("
Distance sum: x = %f y = %f z = %f", $xzy[$key_min][0], $xzy[$key_min][2], $xzy[$key_min][1]); foreach($xzy as $item){ $x[] = $item[0]; $y[] = $item[2]; $z[] = $item[1]; } sort($x); sort($y); sort($z); $cnt = count($x); $half = ($cnt)>>1; $half1 = ($cnt-1)>>1; $med_x = ($x[$half]+$x[$half1])/2; $med_y = ($y[$half]+$y[$half1])/2; $med_z = ($z[$half]+$z[$half1])/2; printf("
Undependent medians: x = %f y = %f z = %f", $med_x, $med_y, $med_z); var_dump($x); var_dump($y); var_dump($z); Результаты: Distance sum: x = -0.763561 y = -3.532885 z = 0.037422 Undependent medians: x = -1.122307 y = -3.571598 z = 0.166906 array (size=10) 0 => float -2.12479996681 1 => float -1.72830736637 2 => float -1.58308649063 3 => float -1.47408938408 4 => float -1.1618667841 5 => float -1.08274650574 6 => float -0.763560533524 7 => float 0.52959895134 8 => float 2.70908093452 9 => float 2.92364525795 array (size=10) 0 => float -4.85960626602 1 => float -4.48901081085 2 => float -4.16520500183 3 => float -3.85271120071 4 => float -3.61031103134 5 => float -3.53288459778 6 => float -3.51698851585 7 => float -1.036495924 8 => float 0.141526937485 9 => float 1.55165171623 array (size=10) 0 => float -1.93620848656 1 => float -0.803761303425 2 => float -0.403136968613 3 => float -0.20462243259 4 => float 0.0374217033386 5 => float 0.296390175819 6 => float 0.345854461193 7 => float 1.11269950867 8 => float 1.23389625549 9 => float 3.04909181595 Оба варианта неплохие. Мне больше нравится медиана. По поводу кластеризации Если ограничиться случаем, когда каждый кластер имеет представителем одну из исходных точек, то множество (набор) точек-представителей кластеров является подмножеством множества исходных точек. Этот набор должен обеспечивать минимум (по всем наборам представителей) из максимумов (по всем исходным точкам) расстояний от этих точек до ближайшего представителя. Но такая задача может быть решена прямым перебором, оптимальность которого очевидна, а проигрыш во времени - не очевиден. Соответствующая программа на PHP: $xzy = [ [ 2.709080934524536 , -1.9362084865570068 , -1.0364959239959717 ], [ 2.923645257949829 , 1.1126995086669922 , 1.5516517162322998 ], [ 0.5295989513397217 , 3.0490918159484863 , 0.1415269374847412 ], [ -1.4740893840789795 , 1.233896255493164 , -3.516988515853882 ], [ -2.124799966812134 , 0.3458544611930847 , -4.8596062660217285 ], [ -1.0827465057373047 , -0.8037613034248352 , -4.489010810852051 ], [ -0.7635605335235596 , 0.03742170333862305 , -3.5328845977783203 ], [ -1.1618667840957642 , -0.4031369686126709 , -3.6103110313415527 ], [ -1.5830864906311035 , 0.296390175819397 , -3.8527112007141113 ], [ -1.7283073663711548 , -0.2046224325895309 , -4.165205001831055 ] ]; function distance_sqr($a, $b){ // квадрат расстояния между двумя точками return ($a[0]-$b[0])*($a[0]-$b[0]) + ($a[1]-$b[1])*($a[1]-$b[1]) + ($a[2]-$b[2])*($a[2]-$b[2]); } function distances($arr){ // матрица расстояний между точками $r = []; foreach($arr as $key1 => $coord1){ $r[$key1] = []; foreach($arr as $key2 => $coord2){ if($key1 <= $key2){ $r[$key1][$key2]= sqrt(distance_sqr($coord1, $coord2)); }else{ $r[$key1][$key2] = $r[$key2][$key1]; } } } return $r; } function combinations($points, $clust){ // генерация наборов представителей $combi_old = [[]]; for($c = 0; $c < $clust; $c++){ $combi = []; foreach($combi_old as $item){ $p0 = empty($item) ? 0 : end($item)+1; for($p = $p0; $p < $points; $p++){ $combi[] = array_merge($item,[$p]); } } $combi_old = $combi; } return $combi; } function clusters($dist, $combi){ // минимакс расстояний до ближайшего кластера $cnt = count($dist); $clusters = [0]; print("
"); foreach($combi as $list){ // наборы кластеров foreach($dist as $key=>$item){ // по каждому элементу массива вычисляются $distances = []; // расстояния до кластеров foreach($list as $index){ $distances[] = $item[$index]; } $d[$key] = min($distances); // расстояние до ближ. кластера } $dd = max($d); // максимум по точкам if(($clusters[0] == 0) || ($dd < $clusters[0])){ // минимум по наборам $clusters = array_merge([$dd], $list); } } return $clusters; } function report($dist, $cluster){ // группировка точек по кластерам $full_clusters = []; foreach($cluster as $clust){ // создаём массивы точек кластера if(is_float($clust)) continue; $full_clusters[$clust] = []; } foreach($dist as $key => $d){ //var_dump($d); $clust_min = -1; $dist_min = null; foreach($full_clusters as $clust => $item){ if(is_null($dist_min) || ($d[$clust] < $dist_min)){ $dist_min = $d[$clust]; $clust_min = $clust; } } $full_clusters[$clust_min][$key] = $dist_min; } //var_dump($full_clusters); return $full_clusters; } print("Координаты точек:
"); foreach($xzy as $c){ printf("[%5.1f,%5.1f, %5.1f] ", $c[0], $c[2], $c[1]); } print("
Матрица расстояний:"); $dist = distances($xzy); foreach($dist as $d){ print("
"); foreach($d as $dd){ printf("%5.1f ", $dd); } } print("
"); $points = count($xzy); for($clust = 1; $clust < 10; $clust++){ $combi = combinations($points, $clust); $cluster = clusters($dist, $combi); print("
Точек: $points Кластеров: $clust. Результат: "); foreach($cluster as $c){ if(is_float($c)) printf("%5.2f ", $c); else printf("% 2d ", $c); } printf("
Нормированная невязка: %5.2f
", $cluster[0]/($points-$clust)); $report = report($dist, $cluster); foreach($report as $key => $full_clust){ print("
Кластер $key: "); foreach($full_clust as $point => $val){ printf("% 2d => % 5.2f ", $point, $val); } } } Результаты: Координаты точек: [ 2.7, -1.0, -1.9] [ 2.9, 1.6, 1.1] [ 0.5, 0.1, 3.0] [ -1.5, -3.5, 1.2] [ -2.1, -4.9, 0.3] [ -1.1, -4.5, -0.8] [ -0.8, -3.5, 0.0] [ -1.2, -3.6, -0.4] [ -1.6, -3.9, 0.3] [ -1.7, -4.2, -0.2] Матрица расстояний: 0.0 4.0 5.6 5.8 6.6 5.3 4.7 4.9 5.6 5.7 4.0 0.0 3.4 6.7 8.2 7.5 6.4 6.8 7.1 7.5 5.6 3.4 0.0 4.5 6.3 6.2 4.9 5.4 5.3 5.9 5.8 6.7 4.5 0.0 1.7 2.3 1.4 1.7 1.0 1.6 6.6 8.2 6.3 1.7 0.0 1.6 1.9 1.7 1.1 1.0 5.3 7.5 6.2 2.3 1.6 0.0 1.3 1.0 1.4 0.9 4.7 6.4 4.9 1.4 1.9 1.3 0.0 0.6 0.9 1.2 4.9 6.8 5.4 1.7 1.7 1.0 0.6 0.0 0.9 0.8 5.6 7.1 5.3 1.0 1.1 1.4 0.9 0.9 0.0 0.6 5.7 7.5 5.9 1.6 1.0 0.9 1.2 0.8 0.6 0.0 Точек: 10 Кластеров: 1. Результат: 6.27 2 Нормированная невязка: 0.70 Кластер 2: 0 => 5.57 1 => 3.39 2 => 0.00 3 => 4.55 4 => 6.27 5 => 6.24 6 => 4.92 7 => 5.37 8 => 5.29 9 => 5.85 Точек: 10 Кластеров: 2. Результат: 4.01 1 3 Нормированная невязка: 0.50 Кластер 1: 0 => 4.01 1 => 0.00 2 => 3.39 Кластер 3: 3 => 0.00 4 => 1.74 5 => 2.29 6 => 1.39 7 => 1.67 8 => 1.00 9 => 1.60 Точек: 10 Кластеров: 3. Результат: 3.39 0 1 3 Нормированная невязка: 0.48 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 2 => 3.39 Кластер 3: 3 => 0.00 4 => 1.74 5 => 2.29 6 => 1.39 7 => 1.67 8 => 1.00 9 => 1.60 Точек: 10 Кластеров: 4. Результат: 1.37 0 1 2 8 Нормированная невязка: 0.23 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 8: 3 => 1.00 4 => 1.14 5 => 1.37 6 => 0.92 7 => 0.85 8 => 0.00 9 => 0.61 Точек: 10 Кластеров: 5. Результат: 1.00 0 1 2 8 9 Нормированная невязка: 0.20 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 8: 3 => 1.00 6 => 0.92 8 => 0.00 Кластер 9: 4 => 0.97 5 => 0.94 7 => 0.82 9 => 0.00 Точек: 10 Кластеров: 6. Результат: 0.97 0 1 2 3 4 7 Нормированная невязка: 0.24 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 3: 3 => 0.00 Кластер 4: 4 => 0.00 Кластер 7: 5 => 0.97 6 => 0.60 7 => 0.00 8 => 0.85 9 => 0.82 Точек: 10 Кластеров: 7. Результат: 0.85 0 1 2 3 4 5 7 Нормированная невязка: 0.28 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 3: 3 => 0.00 Кластер 4: 4 => 0.00 Кластер 5: 5 => 0.00 Кластер 7: 6 => 0.60 7 => 0.00 8 => 0.85 9 => 0.82 Точек: 10 Кластеров: 8. Результат: 0.61 0 1 2 3 4 5 6 8 Нормированная невязка: 0.30 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 3: 3 => 0.00 Кластер 4: 4 => 0.00 Кластер 5: 5 => 0.00 Кластер 6: 6 => 0.00 7 => 0.60 Кластер 8: 8 => 0.00 9 => 0.61 Точек: 10 Кластеров: 9. Результат: 0.60 0 1 2 3 4 5 6 8 9 Нормированная невязка: 0.60 Кластер 0: 0 => 0.00 Кластер 1: 1 => 0.00 Кластер 2: 2 => 0.00 Кластер 3: 3 => 0.00 Кластер 4: 4 => 0.00 Кластер 5: 5 => 0.00 Кластер 6: 6 => 0.00 7 => 0.60 Кластер 8: 8 => 0.00 Кластер 9: 9 => 0.00 Для оценки качества кластеризации минимаксное расстояние делится на разность между количеством точек и количеством кластеров (нормированная невязка). Поскольку расстояния нигде не суммируются, при вычислении матрицы расстояний можно использовать их квадраты. Оценка для случая одного кластера - ещё один вариант ответа на поставленный вопрос.Ответ 2
Привожу в качестве примера алгоритм максимина, который мы когда-то рассматривали в универе, для динамического определения кластеров. Исходные данные: массив точек Алгоритм: Из массива точек выбираются две максимально удаленные друг от друга точки, эти две точки становятся первыми двумя кластерами. Для каждой точки из начального массива высчитывается расстояние до имеющихся кластеров, для каждой точки сохраняется расстояние только до ближайшего кластера (минимальное). Из всех минимальных расстояний выбирается максимальное (максимин) Высчитывается пороговое значение - берется минимальное расстояние между всеми кластерами и делится на два. Если максимин больше чем пороговое значение, то точка становится новым кластером. Повторяем 2-4 пока максимин больше порогового значения, иначе выходим из алгоритма, все кластеры найдены. Для расчета расстояния между точками использовалось расстояние по Евклиду. Точка относится к тому кластеру, расстояние до которого от нее минимально. Получилось адаптировать мою старую наработку под вашу задачу. Реализация не оптимизирована, просто демонстрирует алгоритм, написано на JavaScript. Ссылка на JSFiddle ---Ответ 3
Предисловие Заинтересовала возможность реализации кластеризации в рамках данной задачи. Поскольку ранее этим не занимался, то возможные предложения по улучшению алгоритма приветствуются. Описание алгоритма Алгоритм кластеризации и выделения основной группы объектов основан на построении минимального остовного дерева (minimal spanning tree, MST) для множества несвязных точек в трехмерном пространстве. Построенное дерево является взвешенным, вес ребра - евклидово расстояние между точками в пространстве. На основе полученного дерева определяется средний вес всех рёбер threshold - эта величина используется как пороговое значение принадлежности точки к кластеру (рёбра, чей вес больше этого значения, удаляются). Разбив таким образом MST на отдельные поддеревья (кластеры), определяется кластер с наибольшим количеством объектов, для которого и рассчитывается геометрический центр. Для построения MST используется алгоритм Крускала: Вначале текущее множество рёбер устанавливается пустым. Затем, пока это возможно, проводится следующая операция: из всех рёбер, добавление которых к уже имеющемуся множеству не вызовет появление в нём цикла, выбирается ребро минимального веса и добавляется к уже имеющемуся множеству. Когда таких рёбер больше нет, алгоритм завершён. Для проверки возможной цикличности при соединении рёбер используется union-find data type (система непересекающихся множеств), имеющая следующий интерфейс: /** * Инициализация структуры данных с начальным числом * objectsCount объектов (от 0 до objectsCount - 1). * Каждый объект представляет одно множество. */ public UnionFindStructure(int objectsCount); /** * Возвращает идентификатор множества, в котором расположен объект p. */ public int find(int p); /** * Соединяет множество, содержащее объект p, со множеством, содержащим q. */ public void union(int p, int q); /** * Возвращает количество непересекающихся множеств. */ public int count(); /** * Проверяет расположены ли объекты p и q в одном множестве. */ public boolean connected(int p, int q); Реализация алгоритма Крускала: private static void runKruskalMST(Graph graph) { final Listvertices = graph.getVertices(); final int verticesCount = vertices.size(); // Построение множества всех возможных ребер для графа // (другими словами, соединяем каждую вершину друг с другом). final Set existEdges = new HashSet<>(); final Queue minEdgeQueue = new PriorityQueue<>(graph.getVerticesCount(), new Comparator () { @Override public int compare(Edge o1, Edge o2) { return Double.compare(o1.getWeight(), o2.getWeight()); } }); for (int i = 0; i < verticesCount; ++i) { final Vertex currentVertex = vertices.get(i); for (int j = 0; j < verticesCount; ++j) { if (j == i) { continue; } final Vertex adjacentVertex = vertices.get(j); final long edgeId = Edge.id(currentVertex, adjacentVertex); if (!existEdges.contains(edgeId)) { final double edgeWeight = calcEuclideanDistance(currentVertex.getData(), adjacentVertex.getData()); final Edge edge = new Edge(currentVertex, adjacentVertex, edgeWeight); existEdges.add(edgeId); minEdgeQueue.offer(edge); } } } // Используем union-find data type для того, чтобы // избежать замыкания дерева. final UnionFindStructure uf = new UnionFindStructure(verticesCount); // Выполняем построение MST пока есть доступные связи. // Минимальное остовное дерево считается построенным, если // количество связей равно количеству вершин в графе минус один. while (!minEdgeQueue.isEmpty() && graph.getEdgesCount() < graph.getVerticesCount() - 1) { final Edge e = minEdgeQueue.poll(); final int p = e.getSourceVertex().getOrderNumber(); final int q = e.getTargetVertex().getOrderNumber(); if (!uf.connected(p, q)) { uf.union(p, q); graph.connect(e); } } } Ознакомиться с полным кодом проекта можно на BitBucket. По моим подсчетам асимптотическая сложность алгоритма составляет O(n^2), где n - количество исходных точек (объектов). Результат работы на тестовом наборе данных Количество вершин в графе: 10 Средний вес между вершинами минимального остовного дерева: 1.8751318633294416 === Связи построенного минимального остовного дерева === "P8" => "P7" (0,598945). "P9" => "P10" (0,608075). "P1" => "P2" (4,005045). "P3" => "P4" (4,549114). "P9" => "P4" (1,001752). "P2" => "P3" (3,386669). "P6" => "P10" (0,938385). "P8" => "P10" (0,817417). "P10" => "P5" (0,970785). ======================================================== Разрываем связь: "P1" => "P2" (4,005045). Разрываем связь: "P3" => "P4" (4,549114). Разрываем связь: "P2" => "P3" (3,386669). Количество кластеров: 4 === Кластер #1. Количество объектов: 1 === P1(2,709081, -1,036496, -1,936208) ======================================================== === Кластер #2. Количество объектов: 1 === P2(2,923645, 1,551652, 1,112700) ======================================================== === Кластер #3. Количество объектов: 1 === P3(0,529599, 0,141527, 3,049092) ======================================================== === Кластер #4. Количество объектов: 7 === P9(-1,583086, -3,852711, 0,296390) P6(-1,082747, -4,489011, -0,803761) P8(-1,161867, -3,610311, -0,403137) P4(-1,474089, -3,516989, 1,233896) P7(-0,763561, -3,532885, 0,037422) P10(-1,728307, -4,165205, -0,204622) P5(-2,124800, -4,859606, 0,345854) ======================================================== Основная группа объектов представлена кластером #4 Положение центра основной группы объектов: PivotPoint(-1,416922, -4,003817, 0,071720) Дополнительные источники Intro to Union-Find Data Structure Minimum Spanning Trees, Algorithms, 4th Edition
Комментариев нет:
Отправить комментарий