diff options
Diffstat (limited to 'extern/quadriflow/src/Optimizer.cu')
| -rw-r--r-- | extern/quadriflow/src/Optimizer.cu | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/extern/quadriflow/src/Optimizer.cu b/extern/quadriflow/src/Optimizer.cu new file mode 100644 index 00000000000..8520c3f90a8 --- /dev/null +++ b/extern/quadriflow/src/Optimizer.cu @@ -0,0 +1,281 @@ +#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); +} |