diff options
Diffstat (limited to 'source/blender/editors/uvedit/uvedit_clipboard_graph_iso.cc')
| -rw-r--r-- | source/blender/editors/uvedit/uvedit_clipboard_graph_iso.cc | 419 |
1 files changed, 419 insertions, 0 deletions
diff --git a/source/blender/editors/uvedit/uvedit_clipboard_graph_iso.cc b/source/blender/editors/uvedit/uvedit_clipboard_graph_iso.cc new file mode 100644 index 00000000000..43975d0d382 --- /dev/null +++ b/source/blender/editors/uvedit/uvedit_clipboard_graph_iso.cc @@ -0,0 +1,419 @@ +/* SPDX-License-Identifier: GPL-3.0-or-later + * Copyright (c) 2019 Stefano Quer. + * Additional code, copyright 2022 Blender Foundation. All rights reserved. + * + * Originally 6846114 from https://github.com/stefanoquer/graphISO/blob/master/v3 + * graphISO: Tools to compute the Maximum Common Subgraph between two graphs. + */ + +#include "uvedit_clipboard_graph_iso.hh" + +#include "BLI_assert.h" + +#include "MEM_guardedalloc.h" + +#include <algorithm> +#include <limits.h> + +#define L 0 +#define R 1 +#define LL 2 +#define RL 3 +#define ADJ 4 +#define P 5 +#define W 6 +#define IRL 7 + +#define BDS 8 + +GraphISO::GraphISO(int n) +{ + this->n = n; + label = static_cast<uint *>(MEM_mallocN(n * sizeof *label, __func__)); + adjmat = static_cast<uint8_t **>(MEM_mallocN(n * sizeof *adjmat, __func__)); + + /* \note Allocation of `n * n` bytes total! */ + + for (int i = 0; i < n; i++) { + /* Caution, are you trying to change the representation of adjmat? + * Consider `blender::Vector<std::pair<int, int>> adjmat;` instead. + * Better still is to use a different algorithm. See for example: + * https://www.uni-ulm.de/fileadmin/website_uni_ulm/iui.inst.190/Mitarbeiter/toran/beatcs09.pdf + */ + adjmat[i] = static_cast<unsigned char *>(MEM_callocN(n * sizeof *adjmat[i], __func__)); + } + degree = nullptr; +} + +GraphISO::~GraphISO() +{ + for (int i = 0; i < n; i++) { + MEM_freeN(adjmat[i]); + } + MEM_freeN(adjmat); + MEM_freeN(label); + if (degree) { + MEM_freeN(degree); + } +} + +void GraphISO::add_edge(int v, int w) +{ + BLI_assert(v != w); + adjmat[v][w] = 1; + adjmat[w][v] = 1; +} + +void GraphISO::calculate_degrees() const +{ + if (degree) { + return; + } + degree = static_cast<unsigned int *>(MEM_mallocN(n * sizeof *degree, __func__)); + for (int v = 0; v < n; v++) { + int row_count = 0; + for (int w = 0; w < n; w++) { + if (adjmat[v][w]) { + row_count++; + } + } + degree[v] = row_count; + } +} + +class GraphISO_DegreeCompare { + public: + GraphISO_DegreeCompare(const GraphISO *g) + { + this->g = g; + } + const GraphISO *g; + + bool operator()(int i, int j) const + { + return g->degree[i] < g->degree[j]; + } +}; + +GraphISO *GraphISO::sort_vertices_by_degree() const +{ + calculate_degrees(); + + int *vv = static_cast<int *>(MEM_mallocN(n * sizeof *vv, __func__)); + for (int i = 0; i < n; i++) { + vv[i] = i; + } + /* Currently ordering iso_verts by degree. + * Instead should order iso_verts by frequency of degree. */ + GraphISO_DegreeCompare compare_obj(this); + std::sort(vv, vv + n, compare_obj); + + GraphISO *subg = new GraphISO(n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + subg->adjmat[i][j] = adjmat[vv[i]][vv[j]]; + } + } + for (int i = 0; i < n; i++) { + subg->label[i] = label[vv[i]]; + } + subg->calculate_degrees(); + + MEM_freeN(vv); + return subg; +} + +static void update_incumbent(uint8_t cur[][2], int inc[][2], int cur_pos, int *inc_pos) +{ + if (cur_pos > *inc_pos) { + *inc_pos = cur_pos; + for (int i = 0; i < cur_pos; i++) { + inc[i][L] = cur[i][L]; + inc[i][R] = cur[i][R]; + } + } +} + +static void add_bidomain(uint8_t domains[][BDS], + int *bd_pos, + uint8_t left_i, + uint8_t right_i, + uint8_t left_len, + uint8_t right_len, + uint8_t is_adjacent, + uint8_t cur_pos) +{ + domains[*bd_pos][L] = left_i; + domains[*bd_pos][R] = right_i; + domains[*bd_pos][LL] = left_len; + domains[*bd_pos][RL] = right_len; + domains[*bd_pos][ADJ] = is_adjacent; + domains[*bd_pos][P] = cur_pos; + domains[*bd_pos][W] = UINT8_MAX; + domains[*bd_pos][IRL] = right_len; + (*bd_pos)++; +} + +static int calc_bound(uint8_t domains[][BDS], int bd_pos, int cur_pos) +{ + int bound = 0; + for (int i = bd_pos - 1; i >= 0 && domains[i][P] == cur_pos; i--) { + bound += std::min(domains[i][LL], domains[i][IRL]); + } + return bound; +} + +static int partition(uint8_t *arr, int start, int len, const uint8_t *adjrow) +{ + int i = 0; + for (int j = 0; j < len; j++) { + if (adjrow[arr[start + j]]) { + std::swap(arr[start + i], arr[start + j]); + i++; + } + } + return i; +} + +static void generate_next_domains(uint8_t domains[][BDS], + int *bd_pos, + int cur_pos, + uint8_t *left, + uint8_t *right, + uint8_t v, + uint8_t w, + int inc_pos, + uint8_t **adjmat0, + uint8_t **adjmat1) +{ + int i; + int bd_backup = *bd_pos; + int bound = 0; + uint8_t *bd; + for (i = *bd_pos - 1, bd = &domains[i][L]; i >= 0 && bd[P] == cur_pos - 1; + i--, bd = &domains[i][L]) { + + uint8_t l_len = partition(left, bd[L], bd[LL], adjmat0[v]); + uint8_t r_len = partition(right, bd[R], bd[RL], adjmat1[w]); + + if (bd[LL] - l_len && bd[RL] - r_len) { + add_bidomain(domains, + bd_pos, + bd[L] + l_len, + bd[R] + r_len, + bd[LL] - l_len, + bd[RL] - r_len, + bd[ADJ], + (uint8_t)(cur_pos)); + bound += std::min(bd[LL] - l_len, bd[RL] - r_len); + } + if (l_len && r_len) { + add_bidomain(domains, bd_pos, bd[L], bd[R], l_len, r_len, true, (uint8_t)(cur_pos)); + bound += std::min(l_len, r_len); + } + } + if (cur_pos + bound <= inc_pos) { + *bd_pos = bd_backup; + } +} + +static uint8_t select_next_v(uint8_t *left, uint8_t *bd) +{ + uint8_t min = UINT8_MAX; + uint8_t idx = UINT8_MAX; + if (bd[RL] != bd[IRL]) { + return left[bd[L] + bd[LL]]; + } + for (uint8_t i = 0; i < bd[LL]; i++) { + if (left[bd[L] + i] < min) { + min = left[bd[L] + i]; + idx = i; + } + } + std::swap(left[bd[L] + idx], left[bd[L] + bd[LL] - 1]); + bd[LL]--; + bd[RL]--; + return min; +} + +static uint8_t find_min_value(uint8_t *arr, uint8_t start_idx, uint8_t len) +{ + uint8_t min_v = UINT8_MAX; + for (int i = 0; i < len; i++) { + if (arr[start_idx + i] < min_v) { + min_v = arr[start_idx + i]; + } + } + return min_v; +} + +static void select_bidomain( + uint8_t domains[][BDS], int bd_pos, uint8_t *left, int current_matching_size, bool connected) +{ + int i; + int min_size = INT_MAX; + int min_tie_breaker = INT_MAX; + int best = INT_MAX; + uint8_t *bd; + for (i = bd_pos - 1, bd = &domains[i][L]; i >= 0 && bd[P] == current_matching_size; + i--, bd = &domains[i][L]) { + if (connected && current_matching_size > 0 && !bd[ADJ]) { + continue; + } + int len = bd[LL] > bd[RL] ? bd[LL] : bd[RL]; + if (len < min_size) { + min_size = len; + min_tie_breaker = find_min_value(left, bd[L], bd[LL]); + best = i; + } + else if (len == min_size) { + int tie_breaker = find_min_value(left, bd[L], bd[LL]); + if (tie_breaker < min_tie_breaker) { + min_tie_breaker = tie_breaker; + best = i; + } + } + } + if (best != INT_MAX && best != bd_pos - 1) { + uint8_t tmp[BDS]; + for (i = 0; i < BDS; i++) { + tmp[i] = domains[best][i]; + } + for (i = 0; i < BDS; i++) { + domains[best][i] = domains[bd_pos - 1][i]; + } + for (i = 0; i < BDS; i++) { + domains[bd_pos - 1][i] = tmp[i]; + } + } +} + +static uint8_t select_next_w(uint8_t *right, uint8_t *bd) +{ + uint8_t min = UINT8_MAX; + uint8_t idx = UINT8_MAX; + for (uint8_t i = 0; i < bd[RL] + 1; i++) { + if ((right[bd[R] + i] > bd[W] || bd[W] == UINT8_MAX) && right[bd[R] + i] < min) { + min = right[bd[R] + i]; + idx = i; + } + } + if (idx == UINT8_MAX) { + bd[RL]++; + } + return idx; +} + +static void maximum_common_subgraph_internal( + int incumbent[][2], int *inc_pos, uint8_t **adjmat0, int n0, uint8_t **adjmat1, int n1) +{ + int min = std::min(n0, n1); + + uint8_t cur[min][2]; + uint8_t domains[min * min][BDS]; + uint8_t left[n0], right[n1]; + uint8_t v, w, *bd; + int bd_pos = 0; + for (int i = 0; i < n0; i++) { + left[i] = i; + } + for (int i = 0; i < n1; i++) { + right[i] = i; + } + add_bidomain(domains, &bd_pos, 0, 0, n0, n1, 0, 0); + + int iteration_count = 0; + + while (bd_pos > 0) { + if (iteration_count++ > 10000000) { + /* Unlikely to find a solution, may as well give up. + * Can occur with moderate sized inputs where the graph has lots of symmetry, e.g. a cube + * subdivided 3x times. + */ + return; + } + bd = &domains[bd_pos - 1][L]; + if (calc_bound(domains, bd_pos, bd[P]) + bd[P] <= *inc_pos || + (bd[LL] == 0 && bd[RL] == bd[IRL])) { + bd_pos--; + } + else { + const bool connected = false; + select_bidomain(domains, bd_pos, left, domains[bd_pos - 1][P], connected); + v = select_next_v(left, bd); + if ((bd[W] = select_next_w(right, bd)) != UINT8_MAX) { + w = right[bd[R] + bd[W]]; /* Swap the W after the bottom of the current right domain. */ + right[bd[R] + bd[W]] = right[bd[R] + bd[RL]]; + right[bd[R] + bd[RL]] = w; + bd[W] = w; /* Store the W used for this iteration. */ + cur[bd[P]][L] = v; + cur[bd[P]][R] = w; + update_incumbent(cur, incumbent, bd[P] + (uint8_t)1, inc_pos); + generate_next_domains( + domains, &bd_pos, bd[P] + 1, left, right, v, w, *inc_pos, adjmat0, adjmat1); + } + } + } +} + +static bool check_automorphism(const GraphISO *g0, + const GraphISO *g1, + int solution[][2], + int *solution_length) +{ + if (g0->n != g1->n) { + return false; + } + for (int i = 0; i < g0->n; i++) { + if (g0->label[i] != g1->label[i]) { + return false; + } + for (int j = 0; j < g0->n; j++) { + if (g0->adjmat[i][j] != g1->adjmat[i][j]) { + return false; + } + } + solution[i][0] = i; + solution[i][1] = i; + } + *solution_length = g0->n; + return true; +} + +bool ED_uvedit_clipboard_maximum_common_subgraph(GraphISO *g0_input, + GraphISO *g1_input, + int solution[][2], + int *solution_length) +{ + if (check_automorphism(g0_input, g1_input, solution, solution_length)) { + return true; + } + + int n0 = g0_input->n; + int n1 = g1_input->n; + + int min_size = std::min(n0, n1); + if (min_size >= UINT8_MAX - 2) { + return false; + } + + GraphISO *g0 = g0_input->sort_vertices_by_degree(); + GraphISO *g1 = g1_input->sort_vertices_by_degree(); + + int sol_len = 0; + maximum_common_subgraph_internal(solution, &sol_len, g0->adjmat, n0, g1->adjmat, n1); + *solution_length = sol_len; + + bool result = (sol_len == n0); + if (result) { + for (int i = 0; i < sol_len; i++) { + solution[i][0] = g0->label[solution[i][0]]; /* Index from input. */ + solution[i][1] = g1->label[solution[i][1]]; + } + } + + delete g1; + delete g0; + + return result; +} |