Page MenuHome
Differential D5594 Diff 18116 extern/quadriflow/src/parametrizer-int.cpp

Changeset View

Standalone View

View Options

extern/quadriflow/src/parametrizer-int.cpp

  • This file was added.
#include "parametrizer.hpp"
#include <queue>
#include <unordered_map>
#include <vector>
#include <random>
#include "optimizer.hpp"
namespace qflow {
void Parametrizer::BuildEdgeInfo() {
auto& F = hierarchy.mF;
auto& E2E = hierarchy.mE2E;
edge_diff.clear();
edge_values.clear();
face_edgeIds.resize(F.cols(), Vector3i(-1, -1, -1));
for (int i = 0; i < F.cols(); ++i) {
for (int j = 0; j < 3; ++j) {
int k1 = j, k2 = (j + 1) % 3;
int v1 = F(k1, i);
int v2 = F(k2, i);
DEdge e2(v1, v2);
Vector2i diff2;
int rank2;
if (v1 > v2) {
rank2 = pos_rank(k2, i);
diff2 =
rshift90(Vector2i(-pos_index(k1 * 2, i), -pos_index(k1 * 2 + 1, i)), rank2);
} else {
rank2 = pos_rank(k1, i);
diff2 = rshift90(Vector2i(pos_index(k1 * 2, i), pos_index(k1 * 2 + 1, i)), rank2);
}
int current_eid = i * 3 + k1;
int eid = E2E[current_eid];
int eID1 = face_edgeIds[current_eid / 3][current_eid % 3];
int eID2 = -1;
if (eID1 == -1) {
eID2 = edge_values.size();
edge_values.push_back(e2);
edge_diff.push_back(diff2);
face_edgeIds[i][k1] = eID2;
if (eid != -1) face_edgeIds[eid / 3][eid % 3] = eID2;
} else if (!singularities.count(i)) {
eID2 = face_edgeIds[eid / 3][eid % 3];
edge_diff[eID2] = diff2;
}
}
}
}
void Parametrizer::BuildIntegerConstraints() {
auto& F = hierarchy.mF;
auto& Q = hierarchy.mQ[0];
auto& N = hierarchy.mN[0];
face_edgeOrients.resize(F.cols());
//Random number generator (for shuffling)
std::random_device rd;
std::mt19937 g(rd());
g.seed(hierarchy.rng_seed);
// undirected edge to direct edge
std::vector<std::pair<int, int>> E2D(edge_diff.size(), std::make_pair(-1, -1));
for (int i = 0; i < F.cols(); ++i) {
int v0 = F(0, i);
int v1 = F(1, i);
int v2 = F(2, i);
DEdge e0(v0, v1), e1(v1, v2), e2(v2, v0);
const Vector3i& eid = face_edgeIds[i];
Vector2i variable_id[3];
for (int i = 0; i < 3; ++i) {
variable_id[i] = Vector2i(eid[i] * 2 + 1, eid[i] * 2 + 2);
}
auto index1 =
compat_orientation_extrinsic_index_4(Q.col(v0), N.col(v0), Q.col(v1), N.col(v1));
auto index2 =
compat_orientation_extrinsic_index_4(Q.col(v0), N.col(v0), Q.col(v2), N.col(v2));
int rank1 = (index1.first - index1.second + 4) % 4; // v1 -> v0
int rank2 = (index2.first - index2.second + 4) % 4; // v2 -> v0
int orients[3] = {0}; // == {0, 0, 0}
if (v1 < v0) {
variable_id[0] = -rshift90(variable_id[0], rank1);
orients[0] = (rank1 + 2) % 4;
} else {
orients[0] = 0;
}
if (v2 < v1) {
variable_id[1] = -rshift90(variable_id[1], rank2);
orients[1] = (rank2 + 2) % 4;
} else {
variable_id[1] = rshift90(variable_id[1], rank1);
orients[1] = rank1;
}
if (v2 < v0) {
variable_id[2] = rshift90(variable_id[2], rank2);
orients[2] = rank2;
} else {
variable_id[2] = -variable_id[2];
orients[2] = 2;
}
face_edgeOrients[i] = Vector3i(orients[0], orients[1], orients[2]);
for (int j = 0; j < 3; ++j) {
int eid = face_edgeIds[i][j];
if (E2D[eid].first == -1)
E2D[eid].first = i * 3 + j;
else
E2D[eid].second = i * 3 + j;
}
}
// a face disajoint tree
DisajointOrientTree disajoint_orient_tree = DisajointOrientTree(F.cols());
// merge the whole face graph except for the singularity in which there exists a spanning tree
// which contains consistent orientation
std::vector<int> sharpUE(E2D.size());
for (int i = 0; i < sharp_edges.size(); ++i) {
if (sharp_edges[i]) {
sharpUE[face_edgeIds[i / 3][i % 3]] = 1;
}
}
for (int i = 0; i < E2D.size(); ++i) {
auto& edge_c = E2D[i];
int f0 = edge_c.first / 3;
int f1 = edge_c.second / 3;
if (edge_c.first == -1 || edge_c.second == -1) continue;
if (singularities.count(f0) || singularities.count(f1) || sharpUE[i]) continue;
int orient1 = face_edgeOrients[f0][edge_c.first % 3];
int orient0 = (face_edgeOrients[f1][edge_c.second % 3] + 2) % 4;
disajoint_orient_tree.Merge(f0, f1, orient0, orient1);
}
// merge singularity later
for (auto& f : singularities) {
for (int i = 0; i < 3; ++i) {
if (sharpUE[face_edgeIds[f.first][i]]) continue;
auto& edge_c = E2D[face_edgeIds[f.first][i]];
if (edge_c.first == -1 || edge_c.second == -1) continue;
int v0 = edge_c.first / 3;
int v1 = edge_c.second / 3;
int orient1 = face_edgeOrients[v0][edge_c.first % 3];
int orient0 = (face_edgeOrients[v1][edge_c.second % 3] + 2) % 4;
disajoint_orient_tree.Merge(v0, v1, orient0, orient1);
}
}
for (int i = 0; i < sharpUE.size(); ++i) {
if (sharpUE[i] == 0) continue;
auto& edge_c = E2D[i];
if (edge_c.first == -1 || edge_c.second == -1) continue;
int f0 = edge_c.first / 3;
int f1 = edge_c.second / 3;
int orient1 = face_edgeOrients[f0][edge_c.first % 3];
int orient0 = (face_edgeOrients[f1][edge_c.second % 3] + 2) % 4;
disajoint_orient_tree.Merge(f0, f1, orient0, orient1);
}
// all the face has the same parent. we rotate every face to the space of that parent.
for (int i = 0; i < face_edgeOrients.size(); ++i) {
for (int j = 0; j < 3; ++j) {
face_edgeOrients[i][j] =
(face_edgeOrients[i][j] + disajoint_orient_tree.Orient(i)) % 4;
}
}
std::vector<int> sharp_colors(face_edgeIds.size(), -1);
int num_sharp_component = 0;
// label the connected component connected by non-fixed edges
// we need this because we need sink flow (demand) == source flow (supply) for each component
// rather than global
for (int i = 0; i < sharp_colors.size(); ++i) {
if (sharp_colors[i] != -1) continue;
sharp_colors[i] = num_sharp_component;
std::queue<int> q;
q.push(i);
int counter = 0;
while (!q.empty()) {
int v = q.front();
q.pop();
for (int i = 0; i < 3; ++i) {
int e = face_edgeIds[v][i];
int deid1 = E2D[e].first;
int deid2 = E2D[e].second;
if (deid1 == -1 || deid2 == -1) continue;
if (abs(face_edgeOrients[deid1 / 3][deid1 % 3] -
face_edgeOrients[deid2 / 3][deid2 % 3] + 4) %
4 !=
2 ||
sharpUE[e]) {
continue;
}
for (int k = 0; k < 2; ++k) {
int f = (k == 0) ? E2D[e].first / 3 : E2D[e].second / 3;
if (sharp_colors[f] == -1) {
sharp_colors[f] = num_sharp_component;
q.push(f);
}
}
}
counter += 1;
}
num_sharp_component += 1;
}
{
std::vector<int> total_flows(num_sharp_component);
// check if each component is full-flow
for (int i = 0; i < face_edgeIds.size(); ++i) {
Vector2i diff(0, 0);
for (int j = 0; j < 3; ++j) {
int orient = face_edgeOrients[i][j];
diff += rshift90(edge_diff[face_edgeIds[i][j]], orient);
}
total_flows[sharp_colors[i]] += diff[0] + diff[1];
}
// build "variable"
variables.resize(edge_diff.size() * 2, std::make_pair(Vector2i(-1, -1), 0));
for (int i = 0; i < face_edgeIds.size(); ++i) {
for (int j = 0; j < 3; ++j) {
Vector2i sign = rshift90(Vector2i(1, 1), face_edgeOrients[i][j]);
int eid = face_edgeIds[i][j];
Vector2i index = rshift90(Vector2i(eid * 2, eid * 2 + 1), face_edgeOrients[i][j]);
for (int k = 0; k < 2; ++k) {
auto& p = variables[abs(index[k])];
if (p.first[0] == -1)
p.first[0] = i * 2 + k;
else
p.first[1] = i * 2 + k;
p.second += sign[k];
}
}
}
// fixed variable that might be manually modified.
// modified_variables[component_od][].first = fixed_variable_id
// modified_variables[component_od][].second = 1 if two positive signs -1 if two negative
// signs
std::vector<std::vector<std::pair<int, int>>> modified_variables[2];
for (int i = 0; i < 2; ++i) modified_variables[i].resize(total_flows.size());
for (int i = 0; i < variables.size(); ++i) {
if ((variables[i].first[1] == -1 || variables[i].second != 0) &&
allow_changes[i] == 1) {
int find = sharp_colors[variables[i].first[0] / 2];
int step = std::abs(variables[i].second) % 2;
if (total_flows[find] > 0) {
if (variables[i].second > 0 && edge_diff[i / 2][i % 2] > -1) {
modified_variables[step][find].push_back(std::make_pair(i, -1));
}
if (variables[i].second < 0 && edge_diff[i / 2][i % 2] < 1) {
modified_variables[step][find].push_back(std::make_pair(i, 1));
}
} else if (total_flows[find] < 0) {
if (variables[i].second < 0 && edge_diff[i / 2][i % 2] > -1) {
modified_variables[step][find].push_back(std::make_pair(i, -1));
}
if (variables[i].second > 0 && edge_diff[i / 2][i % 2] < 1) {
modified_variables[step][find].push_back(std::make_pair(i, 1));
}
}
}
}
// uniformly random manually modify variables so that the network has full flow.
for (int i = 0; i < 2; ++i)
for (auto& modified_var : modified_variables[i])
std::shuffle(modified_var.begin(), modified_var.end(), g);
for (int j = 0; j < total_flows.size(); ++j) {
for (int ii = 0; ii < 2; ++ii) {
if (total_flows[j] == 0) continue;
int max_num;
if (ii == 0)
max_num =
std::min(abs(total_flows[j]) / 2, (int)modified_variables[ii][j].size());
else
max_num = std::min(abs(total_flows[j]), (int)modified_variables[ii][j].size());
int dir = (total_flows[j] > 0) ? -1 : 1;
for (int i = 0; i < max_num; ++i) {
auto& info = modified_variables[ii][j][i];
edge_diff[info.first / 2][info.first % 2] += info.second;
if (ii == 0)
total_flows[j] += 2 * dir;
else
total_flows[j] += dir;
}
}
}
}
std::vector<Vector4i> edge_to_constraints(E2D.size() * 2, Vector4i(-1, 0, -1, 0));
for (int i = 0; i < face_edgeIds.size(); ++i) {
for (int j = 0; j < 3; ++j) {
int e = face_edgeIds[i][j];
Vector2i index = rshift90(Vector2i(e * 2 + 1, e * 2 + 2), face_edgeOrients[i][j]);
for (int k = 0; k < 2; ++k) {
int l = abs(index[k]);
int s = index[k] / l;
int ind = l - 1;
int equationID = i * 2 + k;
if (edge_to_constraints[ind][0] == -1) {
edge_to_constraints[ind][0] = equationID;
edge_to_constraints[ind][1] = s;
} else {
edge_to_constraints[ind][2] = equationID;
edge_to_constraints[ind][3] = s;
}
}
}
}
std::vector<std::pair<Vector2i, int>> arcs;
std::vector<int> arc_ids;
DisajointTree tree(face_edgeIds.size() * 2);
for (int i = 0; i < edge_to_constraints.size(); ++i) {
if (allow_changes[i] == 0) continue;
if (edge_to_constraints[i][0] == -1 || edge_to_constraints[i][2] == -1) continue;
if (edge_to_constraints[i][1] == -edge_to_constraints[i][3]) {
int v1 = edge_to_constraints[i][0];
int v2 = edge_to_constraints[i][2];
tree.Merge(v1, v2);
if (edge_to_constraints[i][1] < 0) std::swap(v1, v2);
int current_v = edge_diff[i / 2][i % 2];
arcs.push_back(std::make_pair(Vector2i(v1, v2), current_v));
}
}
tree.BuildCompactParent();
std::vector<int> total_flows(tree.CompactNum());
// check if each component is full-flow
for (int i = 0; i < face_edgeIds.size(); ++i) {
Vector2i diff(0, 0);
for (int j = 0; j < 3; ++j) {
int orient = face_edgeOrients[i][j];
diff += rshift90(edge_diff[face_edgeIds[i][j]], orient);
}
for (int j = 0; j < 2; ++j) {
total_flows[tree.Index(i * 2 + j)] += diff[j];
}
}
// build "variable"
variables.resize(edge_diff.size() * 2);
for (int i = 0; i < variables.size(); ++i) {
variables[i].first = Vector2i(-1, -1);
variables[i].second = 0;
}
for (int i = 0; i < face_edgeIds.size(); ++i) {
for (int j = 0; j < 3; ++j) {
Vector2i sign = rshift90(Vector2i(1, 1), face_edgeOrients[i][j]);
int eid = face_edgeIds[i][j];
Vector2i index = rshift90(Vector2i(eid * 2, eid * 2 + 1), face_edgeOrients[i][j]);
for (int k = 0; k < 2; ++k) {
auto& p = variables[abs(index[k])];
if (p.first[0] == -1)
p.first[0] = i * 2 + k;
else
p.first[1] = i * 2 + k;
p.second += sign[k];
}
}
}
// fixed variable that might be manually modified.
// modified_variables[component_od][].first = fixed_variable_id
// modified_variables[component_od][].second = 1 if two positive signs -1 if two negative signs
std::vector<std::vector<std::pair<int, int>>> modified_variables[2];
for (int i = 0; i < 2; ++i) {
modified_variables[i].resize(total_flows.size());
}
for (int i = 0; i < variables.size(); ++i) {
if ((variables[i].first[1] == -1 || variables[i].second != 0) && allow_changes[i] == 1) {
int find = tree.Index(variables[i].first[0]);
int step = abs(variables[i].second) % 2;
if (total_flows[find] > 0) {
if (variables[i].second > 0 && edge_diff[i / 2][i % 2] > -1) {
modified_variables[step][find].push_back(std::make_pair(i, -1));
}
if (variables[i].second < 0 && edge_diff[i / 2][i % 2] < 1) {
modified_variables[step][find].push_back(std::make_pair(i, 1));
}
} else if (total_flows[find] < 0) {
if (variables[i].second < 0 && edge_diff[i / 2][i % 2] > -1) {
modified_variables[step][find].push_back(std::make_pair(i, -1));
}
if (variables[i].second > 0 && edge_diff[i / 2][i % 2] < 1) {
modified_variables[step][find].push_back(std::make_pair(i, 1));
}
}
}
}
// uniformly random manually modify variables so that the network has full flow.
for (int j = 0; j < 2; ++j) {
for (auto& modified_var : modified_variables[j])
std::shuffle(modified_var.begin(), modified_var.end(), g);
}
for (int j = 0; j < total_flows.size(); ++j) {
for (int ii = 0; ii < 2; ++ii) {
if (total_flows[j] == 0) continue;
int max_num;
if (ii == 0)
max_num = std::min(abs(total_flows[j]) / 2, (int)modified_variables[ii][j].size());
else
max_num = std::min(abs(total_flows[j]), (int)modified_variables[ii][j].size());
int dir = (total_flows[j] > 0) ? -1 : 1;
for (int i = 0; i < max_num; ++i) {
auto& info = modified_variables[ii][j][i];
edge_diff[info.first / 2][info.first % 2] += info.second;
if (ii == 0)
total_flows[j] += 2 * dir;
else
total_flows[j] += dir;
}
}
}
}
void Parametrizer::ComputeMaxFlow() {
hierarchy.DownsampleEdgeGraph(face_edgeOrients, face_edgeIds, edge_diff, allow_changes, 1);
Optimizer::optimize_integer_constraints(hierarchy, singularities, flag_minimum_cost_flow);
hierarchy.UpdateGraphValue(face_edgeOrients, face_edgeIds, edge_diff);
}
} // namespace qflow