~mht/surcut

ref: 16669cc44af1793364ac9aa8db371c55a234ff74 surcut/surcut.cpp -rw-r--r-- 8.9 KiB
16669cc4 — Martin Hafskjold Thoresen Printing and mesh evaluation continues 1 year, 9 months ago
                                                                                
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
#include <igl/per_face_normals.h>

#include "arap-mold.h"
#include "mesh.h"
#include "surcut.h"

std::vector<bool> is_point_above;
std::vector<arap::Moldable::Enum> moldable;
std::vector<arap::Orientation::Enum> orientation;
std::vector<double> dots;
bool regular_arap = false;
float force_scale = 1.0;
bool area_force = false;
float force_dist_cutoff = 1.0;
int iteration_count = 0;
bool stitching = false;
MatrixXd *forces_endpoints_ptr = nullptr;

arap::Data arap_data;
arap::Args arap_args;

Eigen::Vector3d plane_p(0, 0, 0);
Eigen::Vector3d plane_n(0, 1, 0);

// Shoot a `ray` and see if it intersects any of the faces in `mesh`. If it does, return the index
// of the closest vertex to the intersection point. If not, return `-1`.
int ray_intersect(Ray ray, int src_vertex, Mesh &mesh) {
  // XXX(todo): Option to intersect with cut plane. Maybe this isn't a good idea - rays from any
  //            non-moldable face will intersect with the plane!

  Vector3d collision_point;
  double best_dist = std::numeric_limits<double>::max();
  int best_vertex = -1;

  for (int f = 0; f < mesh.faces.rows(); f++) {
    int res = intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.faces,
                                      &collision_point);
    if (res != 1)
      continue;

    double dist_to_collision = std::numeric_limits<double>::max();
    int closest_vertex = -1;

    bool is_own_face = false;
    for (int i = 0; i < 3; i++)
      if (mesh.faces(f, i) == src_vertex)
        is_own_face = true;
    if (is_own_face)
      continue;

    for (int i = 0; i < 3; i++) {
      int v = mesh.faces(f, i);
      double norm = (collision_point - Vector3d(mesh.vertices.row(v))).norm();
      if (norm < dist_to_collision) {
        dist_to_collision = norm;
        closest_vertex = v;
      }
    }

    double dist_to_ray_pos = (Vector3d(mesh.vertices.row(closest_vertex)) - ray.src).norm();
    if (dist_to_ray_pos < best_dist) {
      best_dist = dist_to_ray_pos;
      best_vertex = closest_vertex;
    }
  }
  return best_vertex;
}

// intersect3D_RayTriangle(): find the 3D intersection of a ray with a triangle
//    Input:  a ray R, and a triangle T
//    Output: *I = intersection point (when it exists)
//    Return: -1 = triangle is degenerate (a segment or point)
//             0 =  disjoint (no intersect)
//             1 =  intersect in unique point I1
//             2 =  are in the same plane
// Found at http://geomalgorithms.com/a06-_intersect-2.html#intersect3D_RayTriangle
int intersect3D_RayTriangle(Ray ray, Vector3d triangle_normal, int triangle, MatrixXd &verts,
                            MatrixXi &faces, Vector3d *I) {
#define SMALL_NUM 0.00001
  // get triangle edge vectors and plane normal
  Vector3d v0(verts.row(faces(triangle, 0))), v1(verts.row(faces(triangle, 1))),
      v2(verts.row(faces(triangle, 2)));
  Vector3d u = v1 - v0;
  Vector3d v = v2 - v0;

  Vector3d w0 = ray.src - v0;
  double a = -triangle_normal.dot(w0);
  double b = triangle_normal.dot(ray.dir);
  if (abs(b) < SMALL_NUM) { // ray is  parallel to triangle plane
    if (abs(a) < SMALL_NUM) // ray lies in triangle plane
      return 2;
    else
      return 0; // ray disjoint from plane
  }

  // get intersect point of ray with triangle plane
  double r = a / b;
  if (r < 0.0) // ray goes away from triangle
    return 0;  // => no intersect
  // for a segment, also test if (r > 1.0) => no intersect

  *I = ray.src + r * ray.dir; // intersect point of ray and plane

  // is I inside T?
  double uu, uv, vv, wu, wv, D;
  uu = u.dot(u);
  uv = u.dot(v);
  vv = v.dot(v);
  Vector3d w = *I - v0;
  wu = w.dot(u);
  wv = w.dot(v);
  D = uv * uv - uu * vv;

  // get and test parametric coords
  double s, t;
  s = (uv * wv - vv * wu) / D;
  if (s < 0.0 || s > 1.0) // I is outside T
    return 0;
  t = (uv * wu - uu * wv) / D;
  if (t < 0.0 || (s + t) > 1.0) // I is outside T
    return 0;

  return 1; // I is in T
}

// Perform one ARAP step. Return `true` if the difference between the previous and current state is
// negligible.
StepReturn::Enum step_arap(Mesh &mesh, double scale) {
  iteration_count++;
  Eigen::MatrixXd bc(mesh.boundary.rows(), 3);
  for (int i = 0; i < mesh.boundary.rows(); i++)
    bc.row(i) = mesh.rest_state.row(mesh.boundary(i));

  arap_args.skip_mold_rotations = regular_arap;
  Eigen::MatrixXd prev_vertices = mesh.vertices;
  arap::iter(mesh, arap_args, bc, arap_data, mesh.vertices);
  if (mesh.vertices != mesh.vertices) { /* Contains NaN */
    mesh.vertices = prev_vertices;
    return StepReturn::WouldNaN;
  }
  compute_things(mesh);

  double max_diff = (mesh.vertices - prev_vertices).maxCoeff();
  if (max_diff < scale * 5e-4)
    return StepReturn::Converged;
  return StepReturn::NotDone;
}

// Compute all the things that we need each time the vertex set changes. Here's what we compute:
//  - `TN` (normals)
//  - `is_point_above`
//  - `colors`
void compute_things(Mesh &mesh, double scale, bool only_plane) {
  if (!only_plane) { /* `normals`: depends on vertices */
    igl::per_face_normals(mesh.vertices, mesh.faces, mesh.normals);
  }

  { /* `is_point_above`: depends on vertices and the plane */
    is_point_above.resize(mesh.vertices.rows());
    for (int v_i = 0; v_i < mesh.vertices.rows(); v_i++) {
      Vector3d v(mesh.vertices.row(v_i));
      is_point_above[v_i] = (v - plane_p).dot(plane_n) >= 0.0;
    }
  }

  { /* `moldable`, `orientation`, and `dots`: depends on vertices and the plane
     */
    moldable.resize(mesh.faces.rows());
    orientation.resize(mesh.faces.rows());
    dots.resize(mesh.faces.rows());
    for (int f = 0; f < mesh.faces.rows(); f++) {
      double dot;
      Eigen::Vector3d removal_dir = plane_n;
      moldable[f] =
          arap::is_moldable(removal_dir, is_point_above, mesh.faces, f, mesh.normals.row(f), &dot);
      if (moldable[f] == arap::Moldable::Crossing) {
        orientation[f] = arap::Orientation::Crossing;
      } else {
        if (removal_dir == plane_n) {
          orientation[f] = arap::Orientation::Above;
        } else {
          orientation[f] = arap::Orientation::Below;
        }
      }
      dots[f] = dot;
    }
  }

  { /* `colors`: depends on dot */
    mesh.face_colors.resize(mesh.faces.rows(), 3);
    for (int f = 0; f < mesh.faces.rows(); f++) {
      switch (moldable[f]) {
      case arap::Moldable::Yes: {
        mesh.face_colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
      } break;
      case arap::Moldable::Crossing: {
        mesh.face_colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
      } break;
      default: {
        double r = -dots[f];
        r = sqrt(r);
        mesh.face_colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
      } break;
      }
    }
  }

  if (forces_endpoints_ptr) {
    *forces_endpoints_ptr = mesh.vertices;
  }
  if (stitching) {
    mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), 3);
    std::vector<bool> seen(mesh.vertices.rows());
    for (int f = 0; f < mesh.faces.rows(); f++) {
      auto m = moldable[f];
      // XXX(todo): what do we do about `Crossing`? For now, skip as well.
      if (m == arap::Moldable::Yes)
        continue;
      if (m == arap::Moldable::Crossing) {
        for (int i = 0; i < 3; i++) {
          int v = mesh.faces(f, i);
          seen[v] = true;
          mesh.forces.row(v) = Vector3d::Zero();
        }
        continue;
      }
      Eigen::Vector3d ray_dir =
          (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);

      double area = 1.0;
      if (area_force) {
        Vector3d a = mesh.vertices.row(mesh.faces(f, 1)) - mesh.vertices.row(mesh.faces(f, 0));
        Vector3d b = mesh.vertices.row(mesh.faces(f, 2)) - mesh.vertices.row(mesh.faces(f, 0));
        area = (a.cross(b).norm()) / 2.0;
      }

      for (int i = 0; i < 3; i++) {
        int v = mesh.faces(f, i);
        // Only do a vertex one time; this is okay since the direction is +- plane_n, and not the
        // triangle normal.
        if (seen[v])
          continue;
        seen[v] = true;

        Vector3d vertex = mesh.vertices.row(v);
        Ray ray;
        ray.src = vertex;
        ray.dir = ray_dir;

        int other_i = ray_intersect(ray, v, mesh);
        if (other_i == -1)
          continue;

        Vector3d other_vertex = mesh.vertices.row(other_i);
        auto force_direction = other_vertex - vertex;
        if (force_direction.norm() > scale * force_dist_cutoff)
          continue;

        // XXX(todo): what should the magnitude of the force be? Should probably take a closer look
        // at the dynamics implementation in `arap-mold.cpp`.

        // Symmetric forces - think springs
        Vector3d force =
            force_direction * area * (double)force_scale * force_direction.squaredNorm() / scale;

        mesh.forces.row(v) += force;
        if (forces_endpoints_ptr)
          forces_endpoints_ptr->row(v) = other_vertex;
        mesh.forces.row(other_i) -= force; // Should this depend on area of `other_i`?
      }
    }

    arap_data.f_ext = mesh.forces;
  }
}