Page MenuHome
Differential D5594 Diff 18116 extern/quadriflow/src/Optimizer.cu

Changeset View

Standalone View

View Options

extern/quadriflow/src/Optimizer.cu

  • This file was added.
#include <glm/glm.hpp>
#include <cuda_runtime.h>
#include "AdjacentMatrix.h"
__device__ __host__ glm::dvec3
middle_point(const glm::dvec3 &p0, const glm::dvec3 &n0, const glm::dvec3 &p1, const glm::dvec3 &n1) {
/* How was this derived?
*
* Minimize \|x-p0\|^2 + \|x-p1\|^2, where
* dot(n0, x) == dot(n0, p0)
* dot(n1, x) == dot(n1, p1)
*
* -> Lagrange multipliers, set derivative = 0
* Use first 3 equalities to write x in terms of
* lambda_1 and lambda_2. Substitute that into the last
* two equations and solve for the lambdas. Finally,
* add a small epsilon term to avoid issues when n1=n2.
*/
double n0p0 = glm::dot(n0, p0), n0p1 = glm::dot(n0, p1),
n1p0 = glm::dot(n1, p0), n1p1 = glm::dot(n1, p1),
n0n1 = glm::dot(n0, n1),
denom = 1.0f / (1.0f - n0n1*n0n1 + 1e-4f),
lambda_0 = 2.0f*(n0p1 - n0p0 - n0n1*(n1p0 - n1p1))*denom,
lambda_1 = 2.0f*(n1p0 - n1p1 - n0n1*(n0p1 - n0p0))*denom;
return 0.5 * (p0 + p1) - 0.25 * (n0 * lambda_0 + n1 * lambda_1);
}
__device__ __host__ glm::dvec3
position_round_4(const glm::dvec3 &o, const glm::dvec3 &q,
const glm::dvec3 &n, const glm::dvec3 &p,
double scale) {
double inv_scale = 1.0 / scale;
glm::dvec3 t = glm::cross(n, q);
glm::dvec3 d = p - o;
return o +
q * std::round(glm::dot(q, d) * inv_scale) * scale +
t * std::round(glm::dot(t, d) * inv_scale) * scale;
}
__device__ __host__ glm::dvec3
position_floor_4(const glm::dvec3 &o, const glm::dvec3 &q,
const glm::dvec3 &n, const glm::dvec3 &p,
double scale) {
double inv_scale = 1.0 / scale;
glm::dvec3 t = glm::cross(n,q);
glm::dvec3 d = p - o;
return o +
q * std::floor(glm::dot(q, d) * inv_scale) * scale +
t * std::floor(glm::dot(t, d) * inv_scale) * scale;
}
__device__ __host__ double cudaSignum(double value) {
return std::copysign((double)1, value);
}
__device__ __host__ void
compat_orientation_extrinsic_4(const glm::dvec3 &q0, const glm::dvec3 &n0,
const glm::dvec3 &q1, const glm::dvec3 &n1, glm::dvec3& value1, glm::dvec3& value2) {
const glm::dvec3 A[2] = { q0, glm::cross(n0, q0) };
const glm::dvec3 B[2] = { q1, glm::cross(n1, q1) };
double best_score = -1e10;
int best_a = 0, best_b = 0;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
double score = std::abs(glm::dot(A[i], B[j]));
if (score > best_score + 1e-6) {
best_a = i;
best_b = j;
best_score = score;
}
}
}
const double dp = glm::dot(A[best_a], B[best_b]);
value1 = A[best_a];
value2 = B[best_b] * cudaSignum(dp);
}
__device__ __host__ void
compat_position_extrinsic_4(
const glm::dvec3 &p0, const glm::dvec3 &n0, const glm::dvec3 &q0, const glm::dvec3 &o0,
const glm::dvec3 &p1, const glm::dvec3 &n1, const glm::dvec3 &q1, const glm::dvec3 &o1,
double scale, glm::dvec3& v1, glm::dvec3& v2) {
glm::dvec3 t0 = glm::cross(n0, q0), t1 = glm::cross(n1, q1);
glm::dvec3 middle = middle_point(p0, n0, p1, n1);
glm::dvec3 o0p = position_floor_4(o0, q0, n0, middle, scale);
glm::dvec3 o1p = position_floor_4(o1, q1, n1, middle, scale);
double best_cost = 1e10;
int best_i = -1, best_j = -1;
for (int i = 0; i<4; ++i) {
glm::dvec3 o0t = o0p + (q0 * ((i & 1) * scale) + t0 * (((i & 2) >> 1) * scale));
for (int j = 0; j<4; ++j) {
glm::dvec3 o1t = o1p + (q1 * ((j & 1) * scale) + t1 * (((j & 2) >> 1) * scale));
glm::dvec3 t = o0t - o1t;
double cost = glm::dot(t, t);
if (cost < best_cost) {
best_i = i;
best_j = j;
best_cost = cost;
}
}
}
v1 = o0p + (q0 * ((best_i & 1) * scale) + t0 * (((best_i & 2) >> 1) * scale)),
v2 = o1p + (q1 * ((best_j & 1) * scale) + t1 * (((best_j & 2) >> 1) * scale));
}
__global__
void cudaUpdateOrientation(int* phase, int num_phases, glm::dvec3* N, glm::dvec3* Q, Link* adj, int* adjOffset, int num_adj) {
int pi = blockIdx.x * blockDim.x + threadIdx.x;
// for (int pi = 0; pi < num_phases; ++pi) {
if (pi >= num_phases)
return;
int i = phase[pi];
glm::dvec3 n_i = N[i];
double weight_sum = 0.0f;
glm::dvec3 sum = Q[i];
for (int l = adjOffset[i]; l < adjOffset[i + 1]; ++l) {
Link link = adj[l];
const int j = link.id;
const double weight = link.weight;
if (weight == 0)
continue;
glm::dvec3 n_j = N[j];
glm::dvec3 q_j = Q[j];
glm::dvec3 value1, value2;
compat_orientation_extrinsic_4(sum, n_i, q_j, n_j, value1, value2);
sum = value1 * weight_sum + value2 * weight;
sum -= n_i*glm::dot(n_i, sum);
weight_sum += weight;
double norm = glm::length(sum);
if (norm > 2.93873587705571876e-39f)
sum /= norm;
}
if (weight_sum > 0) {
Q[i] = sum;
}
// }
}
__global__
void cudaPropagateOrientationUpper(glm::dvec3* srcField, glm::ivec2* toUpper, glm::dvec3* N, glm::dvec3* destField, int num_orientation) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
// for (int i = 0; i < num_orientation; ++i) {
if (i >= num_orientation)
return;
for (int k = 0; k < 2; ++k) {
int dest = toUpper[i][k];
if (dest == -1)
continue;
glm::dvec3 q = srcField[i];
glm::dvec3 n = N[dest];
destField[dest] = q - n * glm::dot(n, q);
}
// }
}
__global__
void cudaPropagateOrientationLower(glm::ivec2* toUpper, glm::dvec3* Q, glm::dvec3* N, glm::dvec3* Q_next, glm::dvec3* N_next, int num_toUpper) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
// for (int i = 0; i < num_toUpper; ++i) {
if (i >= num_toUpper)
return;
glm::ivec2 upper = toUpper[i];
glm::dvec3 q0 = Q[upper[0]];
glm::dvec3 n0 = N[upper[0]];
glm::dvec3 q, q1, n1, value1, value2;
if (upper[1] != -1) {
q1 = Q[upper[1]];
n1 = N[upper[1]];
compat_orientation_extrinsic_4(q0, n0, q1, n1, value1, value2);
q = value1 + value2;
}
else {
q = q0;
}
glm::dvec3 n = N_next[i];
q -= glm::dot(n, q) * n;
double len = q.x * q.x + q.y * q.y + q.z * q.z;
if (len > 2.93873587705571876e-39f)
q /= sqrt(len);
Q_next[i] = q;
// }
}
__global__
void cudaUpdatePosition(int* phase, int num_phases, glm::dvec3* N, glm::dvec3* Q, Link* adj, int* adjOffset, int num_adj, glm::dvec3* V, glm::dvec3* O, double scale) {
int pi = blockIdx.x * blockDim.x + threadIdx.x;
// for (int pi = 0; pi < num_phases; ++pi) {
if (pi >= num_phases)
return;
int i = phase[pi];
glm::dvec3 n_i = N[i], v_i = V[i];
glm::dvec3 q_i = Q[i];
glm::dvec3 sum = O[i];
double weight_sum = 0.0f;
for (int l = adjOffset[i]; l < adjOffset[i + 1]; ++l) {
Link link = adj[l];
int j = link.id;
const double weight = link.weight;
if (weight == 0)
continue;
glm::dvec3 n_j = N[j], v_j = V[j];
glm::dvec3 q_j = Q[j], o_j = O[j];
glm::dvec3 v1, v2;
compat_position_extrinsic_4(
v_i, n_i, q_i, sum, v_j, n_j, q_j, o_j, scale, v1, v2);
sum = v1*weight_sum +v2*weight;
weight_sum += weight;
if (weight_sum > 2.93873587705571876e-39f)
sum /= weight_sum;
sum -= glm::dot(n_i, sum - v_i)*n_i;
}
if (weight_sum > 0) {
O[i] = position_round_4(sum, q_i, n_i, v_i, scale);
}
// }
}
__global__
void cudaPropagatePositionUpper(glm::dvec3* srcField, glm::ivec2* toUpper, glm::dvec3* N, glm::dvec3* V, glm::dvec3* destField, int num_position) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
// for (int i = 0; i < num_position; ++i) {
if (i >= num_position)
return;
for (int k = 0; k < 2; ++k) {
int dest = toUpper[i][k];
if (dest == -1)
continue;
glm::dvec3 o = srcField[i], n = N[dest], v = V[dest];
o -= n * glm::dot(n, o - v);
destField[dest] = o;
}
// }
}
void UpdateOrientation(int* phase, int num_phases, glm::dvec3* N, glm::dvec3* Q, Link* adj, int* adjOffset, int num_adj) {
cudaUpdateOrientation << <(num_phases + 255) / 256, 256 >> >(phase, num_phases, N, Q, adj, adjOffset, num_adj);
// cudaUpdateOrientation(phase, num_phases, N, Q, adj, adjOffset, num_adj);
}
void PropagateOrientationUpper(glm::dvec3* srcField, int num_orientation, glm::ivec2* toUpper, glm::dvec3* N, glm::dvec3* destField) {
cudaPropagateOrientationUpper << <(num_orientation + 255) / 256, 256 >> >(srcField, toUpper, N, destField, num_orientation);
// cudaPropagateOrientationUpper(srcField, toUpper, N, destField, num_orientation);
}
void PropagateOrientationLower(glm::ivec2* toUpper, glm::dvec3* Q, glm::dvec3* N, glm::dvec3* Q_next, glm::dvec3* N_next, int num_toUpper) {
cudaPropagateOrientationLower << <(num_toUpper + 255) / 256, 256 >> >(toUpper, Q, N, Q_next, N_next, num_toUpper);
// cudaPropagateOrientationLower(toUpper, Q, N, Q_next, N_next, num_toUpper);
}
void UpdatePosition(int* phase, int num_phases, glm::dvec3* N, glm::dvec3* Q, Link* adj, int* adjOffset, int num_adj, glm::dvec3* V, glm::dvec3* O, double scale) {
cudaUpdatePosition << <(num_phases + 255) / 256, 256 >> >(phase, num_phases, N, Q, adj, adjOffset, num_adj, V, O, scale);
// cudaUpdatePosition(phase, num_phases, N, Q, adj, adjOffset, num_adj, V, O, scale);
}
void PropagatePositionUpper(glm::dvec3* srcField, int num_position, glm::ivec2* toUpper, glm::dvec3* N, glm::dvec3* V, glm::dvec3* destField) {
cudaPropagatePositionUpper << <(num_position + 255) / 256, 256 >> >(srcField, toUpper, N, V, destField, num_position);
// cudaPropagatePositionUpper(srcField, toUpper, N, V, destField, num_position);
}