voronoi3d_algorithm.h 76.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

20
21
#ifndef SWIFT_VORONOIXD_ALGORITHM_H
#define SWIFT_VORONOIXD_ALGORITHM_H
22

23
#include <float.h>
24
#include <math.h>
25
26
#include <stdio.h>
#include <stdlib.h>
27
#include <string.h>
28
#include "error.h"
29
30
31
#include "inline.h"
#include "voronoi3d_cell.h"

32
/* For debugging purposes */
33
//#define LOOP_CHECK 1000
34

35
36
37
38
39
40
#ifdef LOOP_CHECK
/* We need to do the trickery below to get a unique counter for each call to the
   macro. This only works if the macro is never called twice on the same line.
 */
#define MERGE(a, b) a##b
#define LOOPCOUNTER_NAME(line) MERGE(loopcount, line)
41
42

/**
43
 * @brief Increase the given counter variable and check if it is still valid.
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
 *
 * @param counter Counter to increase.
 * @param line_number Line number where the while is called.
 * @return 1 if the counter is still valid, 0 otherwise.
 */
__attribute__((always_inline)) INLINE int check_counter(int *counter,
                                                        int line_number) {
  ++(*counter);
  if ((*counter) == LOOP_CHECK) {
    error("Number of iterations reached maximum (=%i) in while on line %i!",
          LOOP_CHECK, line_number);
  }
  return 1;
}

59
60
61
62
63
/* safewhile is a wrapper around a while that adds a unique counter variable to
   the loop that is increased by 1 for each time the loop is executed, and
   causes the code to crash if this number exceeds a given value.
   We use this to quickly enable or disable number of iterations checks for a
   large number of while loops */
64
65
66
#define safewhile(condition)          \
  int LOOPCOUNTER_NAME(__LINE__) = 0; \
  while (check_counter(&LOOPCOUNTER_NAME(__LINE__), __LINE__) && (condition))
67

68
#else /* LOOP_CHECK */
69
70

/* If LOOP_CHECK is not defined, safewhile and while are EXACTLY the same */
71
#define safewhile(condition) while (condition)
72

73
74
75
76
77
#endif /* LOOP_CHECK */

/* This flag activates a number of expensive geometrical checks that help
   finding bugs. */
//#define VORONOI3D_EXPENSIVE_CHECKS
78

79
80
/* Tolerance parameter used to decide when to use more precise geometric
   criteria */
81
#define VORONOI3D_TOLERANCE 1.e-6f
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

/* Box boundary flags used to signal cells neighbouring the box boundary
   These values correspond to the top range of possible 64-bit integers, and
   we make the strong assumption that there will never be a particle that has
   one of these as particle ID. */
#define VORONOI3D_BOX_FRONT 18446744073709551600llu
#define VORONOI3D_BOX_BACK 18446744073709551601llu
#define VORONOI3D_BOX_TOP 18446744073709551602llu
#define VORONOI3D_BOX_BOTTOM 18446744073709551603llu
#define VORONOI3D_BOX_LEFT 18446744073709551604llu
#define VORONOI3D_BOX_RIGHT 18446744073709551605llu

/*******************************************************************************
 * 3D specific methods
 *
 * Most of these methods are based on the source code of voro++:
 *  http://math.lbl.gov/voro++/
 ******************************************************************************/

101
102
103
104
105
106
107
108
109
/**
 * @brief Print the given cell to the stderr in a format that can be easily
 * plotted using gnuplot.
 *
 * This method prints to the stderr instead of stdout to make it possible to use
 * it right before crashing the code.
 *
 * @param c Voronoi cell to print.
 */
110
__attribute__((always_inline)) INLINE void voronoi_print_gnuplot_c(
111
    const struct voronoi_cell *c) {
112
113

  int i, j, v;
114
  const double *x = c->x;
115
116
117

  fprintf(stderr, "%g\t%g\t%g\n\n", x[0], x[1], x[2]);

118
119
  for (i = 0; i < c->nvert; ++i) {
    for (j = 0; j < c->orders[i]; ++j) {
120
      v = c->edges[c->offsets[i] + j];
121
122
123
      if (v < 0) {
        v = -v - 1;
      }
124
125
126
127
128
129
130
131
132
      fprintf(stderr, "%g\t%g\t%g\n", c->vertices[3 * i + 0] + x[0],
              c->vertices[3 * i + 1] + x[1], c->vertices[3 * i + 2] + x[2]);
      fprintf(stderr, "%g\t%g\t%g\n\n", c->vertices[3 * v + 0] + x[0],
              c->vertices[3 * v + 1] + x[1], c->vertices[3 * v + 2] + x[2]);
    }
  }
  fprintf(stderr, "\n");
}

133
134
135
136
137
138
/**
 * @brief Print the contents of a 3D Voronoi cell
 *
 * @param cell 3D Voronoi cell
 */
__attribute__((always_inline)) INLINE void voronoi_print_cell(
139
    const struct voronoi_cell *cell) {
140
141
142

  int i, j;

143
144
  fprintf(stderr, "x: %g %g %g\n", cell->x[0], cell->x[1], cell->x[2]);
  fprintf(stderr, "nvert: %i\n", cell->nvert);
145

146
  for (i = 0; i < cell->nvert; ++i) {
147
    fprintf(stderr, "%i: %g %g %g (%i)\n", i, cell->vertices[3 * i],
148
149
            cell->vertices[3 * i + 1], cell->vertices[3 * i + 2],
            cell->orders[i]);
150
    for (j = 0; j < cell->orders[i]; ++j) {
151
      fprintf(stderr, "%i (%i)\n", cell->edges[cell->offsets[i] + j],
152
153
154
              cell->edgeindices[cell->offsets[i] + j]);
    }
  }
155
  fprintf(stderr, "\n");
156
157
158
159
}

/**
 * @brief Get the index of the vertex pointed to by the given edge of the given
160
161
162
163
164
165
 * vertex.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @return Index of the vertex on the other side of the edge.
166
167
 */
__attribute__((always_inline)) INLINE int voronoi_get_edge(
168
    const struct voronoi_cell *c, int vertex, int edge) {
169
170
171
172
  return c->edges[c->offsets[vertex] + edge];
}

/**
173
174
175
176
177
178
179
180
181
182
183
184
 * @brief Get the index of the given edge in the edge list of the vertex on the
 * other side of the edge of the given vertex.
 *
 * Suppose that the given vertex has edges [edge1, edge2, given_edge], and that
 * the vertex on the other side of given_edge has edges [edge1, given_edge,
 * edge2], then this method returns 1.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @return Index of that edge in the edge list of the vertex on the other side
 * of the edge.
185
186
 */
__attribute__((always_inline)) INLINE int voronoi_get_edgeindex(
187
    const struct voronoi_cell *c, int vertex, int edge) {
188
189
190
191
  return c->edgeindices[c->offsets[vertex] + edge];
}

/**
192
193
194
195
196
197
198
 * @brief Set the index of the vertex on the other side of the given edge of the
 * given vertex.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of the cell.
 * @param edge Edge of that vertex.
 * @param value Index of the vertex on the other side of that edge.
199
200
201
202
203
204
205
 */
__attribute__((always_inline)) INLINE void voronoi_set_edge(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->edges[c->offsets[vertex] + edge] = value;
}

/**
206
207
208
209
210
211
212
213
214
215
216
217
218
 * @brief Set the index of the given edge in the edge list of the vertex on the
 * other side of the edge of the given vertex.
 *
 * Suppose that the given vertex has edges [edge1, edge2, given_edge], and we
 * want to tell this method that the vertex on the other side of given_edge has
 * edges [edge1, given_edge, edge2], then we need to pass on a value of 1 to
 * this method.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @param value Index of that edge in the edge list of the vertex on the other
 * side of the edge.
219
220
221
222
223
224
225
 */
__attribute__((always_inline)) INLINE void voronoi_set_edgeindex(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->edgeindices[c->offsets[vertex] + edge] = value;
}

/**
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
 * @brief Get the neighbour for the given edge of the given vertex.
 *
 * An edge is shared by two faces, and each face has a neighbour. Luckily, each
 * edge also has two endpoints, so we can get away with storing only one
 * neighbour per endpoint of an edge. We have complete freedom in choosing which
 * neighbour to store in which endpoint, but we need to be consistent about it.
 * Here we use the following convention: if we take a vector pointing away from
 * the given vertex along the given edge direction, then we store the neighbour
 * that corresponds to the face to the right if looking to the cell from the
 * outside. This is the face that you encounter first when rotating counter-
 * clockwise around that vector, starting from outside the cell.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @return Index of the neighbour corresponding to that edge and vertex.
242
243
 */
__attribute__((always_inline)) INLINE int voronoi_get_ngb(
244
    const struct voronoi_cell *c, int vertex, int edge) {
245
246
247
248
  return c->ngbs[c->offsets[vertex] + edge];
}

/**
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
 * @brief Set the neighbour for the given edge of the given vertex.
 *
 * An edge is shared by two faces, and each face has a neighbour. Luckily, each
 * edge also has two endpoints, so we can get away with storing only one
 * neighbour per endpoint of an edge. We have complete freedom in choosing which
 * neighbour to store in which endpoint, but we need to be consistent about it.
 * Here we use the following convention: if we take a vector pointing away from
 * the given vertex along the given edge direction, then we store the neighbour
 * that corresponds to the face to the right if looking to the cell from the
 * outside. This is the face that you encounter first when rotating counter-
 * clockwise around that vector, starting from outside the cell.
 *
 * @param c 3D Voronoi cell.
 * @param vertex Index of a vertex of that cell.
 * @param edge Edge of that vertex.
 * @param value Index of the neighbour corresponding to that edge and vertex.
265
266
267
268
269
270
 */
__attribute__((always_inline)) INLINE void voronoi_set_ngb(
    struct voronoi_cell *c, int vertex, int edge, int value) {
  c->ngbs[c->offsets[vertex] + edge] = value;
}

271
/**
272
273
274
275
276
 * @brief Check if the 3D Voronoi cell is still consistent.
 *
 * A cell is consistent if its edges are consistent, i.e. if edge e of vertex v1
 * points to vertex v2, then v2 should have an edge that points to v1 as well,
 * and then the edge index of vertex v1 should contain the index of that edge
277
278
279
 * in the edge list of v2. We also check if all vertices have orders of at least
 * 3, and if all vertices are actually part of the vertex list.
 * Oh, and we check if the cell actually has vertices.
280
281
282
283
 *
 * @param cell 3D Voronoi cell to check
 */
__attribute__((always_inline)) INLINE void voronoi_check_cell_consistency(
284
    const struct voronoi_cell *c) {
285
286
287

  int i, j, e, l, m;

288
289
290
291
  if (c->nvert < 4) {
    error("Found cell with only %i vertices!", c->nvert);
  }

292
  for (i = 0; i < c->nvert; ++i) {
293
294
295
296
    if (c->orders[i] < 3) {
      voronoi_print_cell(c);
      error("Found cell with vertex of order %i!", c->orders[i]);
    }
297
    for (j = 0; j < c->orders[i]; ++j) {
298
      e = voronoi_get_edge(c, i, j);
299
300
301
302
      if (e >= c->nvert) {
        voronoi_print_cell(c);
        error("Found cell with edges that lead to non-existing vertices!");
      }
303
304
305
      if (e < 0) {
        continue;
      }
306
307
308
      l = voronoi_get_edgeindex(c, i, j);
      m = voronoi_get_edge(c, e, l);
      if (m != i) {
309
        /* voronoi_print_gnuplot_c(c); */
310
        voronoi_print_cell(c);
311
        fprintf(stderr, "i: %i, j: %i, e: %i, l: %i, m: %i\n", i, j, e, l, m);
312
313
314
315
316
317
        error("Cell inconsistency!");
      }
    }
  }
}

318
319
/**
 * @brief Check if the given vertex is above, below or on the cutting plane
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
 * defined by the given parameters.
 *
 * @param v Coordinates of a cell vertex, relative w.r.t. the position of the
 * generator of the cell.
 * @param dx Half of the relative distance vector between the position of the
 * generator of the cell and the position of the neighbouring cell that
 * generates the cutting plane, pointing from the generator position to the
 * cutting plane.
 * @param r2 Squared length of dx.
 * @param test Variable to store the result of the geometric test in, which
 * corresponds to the projected distance between the generator and the vertex
 * along dx.
 * @param teststack Stack to store the results of the N last tests in (for
 * debugging purposes only).
 * @param teststack_size Next available field in the teststack, is reset to 0 if
 * the teststack is full (so the N+1th results is overwritten; for debugging
 * purposes only).
 * @return Result of the test: -1 if the vertex is below the cutting plane, +1
 * if it is above, and 0 if it is on the cutting plane.
339
340
 */
__attribute__((always_inline)) INLINE int voronoi_test_vertex(
341
    const float *v, const float *dx, float r2, float *test, float *teststack,
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    int *teststack_size) {

  *test = v[0] * dx[0] + v[1] * dx[1] + v[2] * dx[2] - r2;

  teststack[*teststack_size] = *test;
  *teststack_size = *teststack_size + 1;
  if (*teststack_size == 2 * VORONOI3D_MAXNUMVERT) {
    *teststack_size = 0;
  }

  if (*test < -VORONOI3D_TOLERANCE) {
    return -1;
  }
  if (*test > VORONOI3D_TOLERANCE) {
    return 1;
  }
  return 0;
}

/**
362
363
364
 * @brief Initialize the cell as a cube that spans the entire simulation box.
 *
 * @param c 3D Voronoi cell to initialize.
365
366
 * @param anchor Anchor of the simulation box.
 * @param side Side lengths of the simulation box.
367
368
 */
__attribute__((always_inline)) INLINE void voronoi_initialize(
369
    struct voronoi_cell *cell, const double *anchor, const double *side) {
370
371
372
373

  cell->nvert = 8;

  /* (0, 0, 0) -- 0 */
374
375
376
  cell->vertices[0] = anchor[0] - cell->x[0];
  cell->vertices[1] = anchor[1] - cell->x[1];
  cell->vertices[2] = anchor[2] - cell->x[2];
377
378

  /* (0, 0, 1)-- 1 */
379
380
381
  cell->vertices[3] = anchor[0] - cell->x[0];
  cell->vertices[4] = anchor[1] - cell->x[1];
  cell->vertices[5] = anchor[2] + side[2] - cell->x[2];
382
383

  /* (0, 1, 0) -- 2 */
384
385
386
  cell->vertices[6] = anchor[0] - cell->x[0];
  cell->vertices[7] = anchor[1] + side[1] - cell->x[1];
  cell->vertices[8] = anchor[2] - cell->x[2];
387
388

  /* (0, 1, 1) -- 3 */
389
390
391
  cell->vertices[9] = anchor[0] - cell->x[0];
  cell->vertices[10] = anchor[1] + side[1] - cell->x[1];
  cell->vertices[11] = anchor[2] + side[2] - cell->x[2];
392
393

  /* (1, 0, 0) -- 4 */
394
395
396
  cell->vertices[12] = anchor[0] + side[0] - cell->x[0];
  cell->vertices[13] = anchor[1] - cell->x[1];
  cell->vertices[14] = anchor[2] - cell->x[2];
397
398

  /* (1, 0, 1) -- 5 */
399
400
401
  cell->vertices[15] = anchor[0] + side[0] - cell->x[0];
  cell->vertices[16] = anchor[1] - cell->x[1];
  cell->vertices[17] = anchor[2] + side[2] - cell->x[2];
402
403

  /* (1, 1, 0) -- 6 */
404
405
406
  cell->vertices[18] = anchor[0] + side[0] - cell->x[0];
  cell->vertices[19] = anchor[1] + side[1] - cell->x[1];
  cell->vertices[20] = anchor[2] - cell->x[2];
407
408

  /* (1, 1, 1) -- 7 */
409
410
411
  cell->vertices[21] = anchor[0] + side[0] - cell->x[0];
  cell->vertices[22] = anchor[1] + side[1] - cell->x[1];
  cell->vertices[23] = anchor[2] + side[2] - cell->x[2];
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532

  cell->orders[0] = 3;
  cell->orders[1] = 3;
  cell->orders[2] = 3;
  cell->orders[3] = 3;
  cell->orders[4] = 3;
  cell->orders[5] = 3;
  cell->orders[6] = 3;
  cell->orders[7] = 3;

  /* edges are ordered counterclockwise w.r.t. a vector pointing from the
     cell generator to the vertex
     (0, 0, 0) corner */
  cell->offsets[0] = 0;
  cell->edges[0] = 1;
  cell->edges[1] = 2;
  cell->edges[2] = 4;
  cell->edgeindices[0] = 0;
  cell->edgeindices[1] = 2;
  cell->edgeindices[2] = 0;

  /* (0, 0, 1) corner */
  cell->offsets[1] = 3;
  cell->edges[3] = 0;
  cell->edges[4] = 5;
  cell->edges[5] = 3;
  cell->edgeindices[3] = 0;
  cell->edgeindices[4] = 2;
  cell->edgeindices[5] = 1;

  /* (0, 1, 0) corner */
  cell->offsets[2] = 6;
  cell->edges[6] = 3;
  cell->edges[7] = 6;
  cell->edges[8] = 0;
  cell->edgeindices[6] = 0;
  cell->edgeindices[7] = 0;
  cell->edgeindices[8] = 1;

  /* (0, 1, 1) corner */
  cell->offsets[3] = 9;
  cell->edges[9] = 2;
  cell->edges[10] = 1;
  cell->edges[11] = 7;
  cell->edgeindices[9] = 0;
  cell->edgeindices[10] = 2;
  cell->edgeindices[11] = 0;

  /* (1, 0, 0) corner */
  cell->offsets[4] = 12;
  cell->edges[12] = 0;
  cell->edges[13] = 6;
  cell->edges[14] = 5;
  cell->edgeindices[12] = 2;
  cell->edgeindices[13] = 2;
  cell->edgeindices[14] = 0;

  /* (1, 0, 1) corner */
  cell->offsets[5] = 15;
  cell->edges[15] = 4;
  cell->edges[16] = 7;
  cell->edges[17] = 1;
  cell->edgeindices[15] = 2;
  cell->edgeindices[16] = 1;
  cell->edgeindices[17] = 1;

  /* (1, 1, 0) corner */
  cell->offsets[6] = 18;
  cell->edges[18] = 2;
  cell->edges[19] = 7;
  cell->edges[20] = 4;
  cell->edgeindices[18] = 1;
  cell->edgeindices[19] = 2;
  cell->edgeindices[20] = 1;

  /* (1, 1, 1) corner */
  cell->offsets[7] = 21;
  cell->edges[21] = 3;
  cell->edges[22] = 5;
  cell->edges[23] = 6;
  cell->edgeindices[21] = 2;
  cell->edgeindices[22] = 1;
  cell->edgeindices[23] = 1;

  /* ngbs[3*i+j] is the neighbour corresponding to the plane clockwise of
     edge j of vertex i (when going from edge j to vertex i)
     we set them to a ridiculously large value to be able to track faces without
     neighbour */
  cell->ngbs[0] = VORONOI3D_BOX_FRONT;  /* (000) - (001) */
  cell->ngbs[1] = VORONOI3D_BOX_LEFT;   /* (000) - (010) */
  cell->ngbs[2] = VORONOI3D_BOX_BOTTOM; /* (000) - (100) */

  cell->ngbs[3] = VORONOI3D_BOX_LEFT;  /* (001) - (000) */
  cell->ngbs[4] = VORONOI3D_BOX_FRONT; /* (001) - (101) */
  cell->ngbs[5] = VORONOI3D_BOX_TOP;   /* (001) - (011) */

  cell->ngbs[6] = VORONOI3D_BOX_LEFT;   /* (010) - (011) */
  cell->ngbs[7] = VORONOI3D_BOX_BACK;   /* (010) - (110) */
  cell->ngbs[8] = VORONOI3D_BOX_BOTTOM; /* (010) - (000) */

  cell->ngbs[9] = VORONOI3D_BOX_BACK;  /* (011) - (010) */
  cell->ngbs[10] = VORONOI3D_BOX_LEFT; /* (011) - (001) */
  cell->ngbs[11] = VORONOI3D_BOX_TOP;  /* (011) - (111) */

  cell->ngbs[12] = VORONOI3D_BOX_FRONT;  /* (100) - (000) */
  cell->ngbs[13] = VORONOI3D_BOX_BOTTOM; /* (100) - (110) */
  cell->ngbs[14] = VORONOI3D_BOX_RIGHT;  /* (100) - (101) */

  cell->ngbs[15] = VORONOI3D_BOX_FRONT; /* (101) - (100) */
  cell->ngbs[16] = VORONOI3D_BOX_RIGHT; /* (101) - (111) */
  cell->ngbs[17] = VORONOI3D_BOX_TOP;   /* (101) - (001) */

  cell->ngbs[18] = VORONOI3D_BOX_BOTTOM; /* (110) - (010) */
  cell->ngbs[19] = VORONOI3D_BOX_BACK;   /* (110) - (111) */
  cell->ngbs[20] = VORONOI3D_BOX_RIGHT;  /* (110) - (100) */

  cell->ngbs[21] = VORONOI3D_BOX_BACK;  /* (111) - (011) */
  cell->ngbs[22] = VORONOI3D_BOX_TOP;   /* (111) - (101) */
  cell->ngbs[23] = VORONOI3D_BOX_RIGHT; /* (111) - (110) */
}

533
/**
534
 * @brief Find an edge of the voronoi_cell that intersects the cutting plane.
535
 *
536
537
538
539
 * There is a large number of possible paths through this method, each of which
 * is covered by a separate unit test in testVoronoi3D. Paths have been numbered
 * in the inline comments to help identify them.
 *
540
 * @param c 3D Voronoi cell.
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
 * @param dx Vector pointing from pj to the midpoint of the line segment between
 * pi and pj.
 * @param r2 Squared length of dx.
 * @param u Projected distance between the plane and the closest vertex above
 * the plane, along dx.
 * @param up Index of the closest vertex above the plane.
 * @param us Index of the edge of vertex up that intersects the plane.
 * @param uw Result of the last test_vertex call for vertex up.
 * @param l Projected distance between the plane and the closest vertex below
 * the plane, along dx.
 * @param lp Index of the closest vertex below the plane.
 * @param ls Index of the edge of vertex lp that intersects the plane.
 * @param lw Result of the last test_vertex call for vertex lp.
 * @param q Projected distance between the plane and a test vertex, along dx.
 * @param qp Index of the test vertex.
 * @param qs Index of the edge of the test vertex that is connected to up.
 * @param qw Result of the last test_vertex call involving qp.
558
559
 * @return A negative value if an error occurred, 0 if the plane does not
 * intersect the cell, 1 if nothing special happened and 2 if we have a
560
 * complicated setup.
561
562
 */
__attribute__((always_inline)) INLINE int voronoi_intersect_find_closest_vertex(
563
564
565
    struct voronoi_cell *c, const float *dx, float r2, float *u, int *up,
    int *us, int *uw, float *l, int *lp, int *ls, int *lw, float *q, int *qp,
    int *qs, int *qw) {
566

567
568
  /* stack to store all vertices that have already been tested (debugging
     only) */
569
  float teststack[2 * VORONOI3D_MAXNUMVERT];
570
  /* size of the used part of the stack */
571
  int teststack_size = 0;
572
  /* flag signalling a complicated setup */
573
574
  int complicated;

575
576
  /* test the first vertex: uw = -1 if it is below the plane, 1 if it is above
     0 if it is very close to the plane, and things become complicated... */
577
578
579
580
581
582
583
584
585
586
587
  *uw = voronoi_test_vertex(&c->vertices[0], dx, r2, u, teststack,
                            &teststack_size);
  *up = 0;
  complicated = 0;
  if ((*uw) == 0) {

    /* PATH 0 */
    complicated = 1;

  } else {

588
    /* two options: either the vertex is above or below the plane */
589
590
591
592
593

    if ((*uw) == 1) {

      /* PATH 1 */

594
595
596
      /* above: try to find a vertex below
         we test all edges of the current vertex stored in up (vertex 0) until
         we either find one below the plane or closer to the plane */
597
598
599
600
      *lp = voronoi_get_edge(c, (*up), 0);
      *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l, teststack,
                                &teststack_size);
      *us = 1;
601
602
603
      /* Not in while: PATH 1.0 */
      /* somewhere in while: PATH 1.1 */
      /* last valid option of while: PATH 1.2 */
604
605
606
607
      safewhile((*us) < c->orders[(*up)] && (*l) >= (*u)) {
        *lp = voronoi_get_edge(c, (*up), (*us));
        *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l, teststack,
                                  &teststack_size);
608
        ++(*us);
609
      }
610
      /* we increased us too much, correct this */
611
      --(*us);
612
613
      if ((*l) >= (*u)) {
        /* PATH 1.3 */
614
615
616
617
618
        /* up is the closest vertex to the plane, but is above the plane
           since the entire cell is convex, up is the closest vertex of all
           vertices of the cell
           this means the entire cell is supposedly above the plane, which is
           impossible */
619
620
621
622
623
624
        message(
            "Cell completely gone! This should not happen. (l >= u, l = %g, u "
            "= %g)",
            (*l), (*u));
        return -1;
      }
625
626
      /* we know that lp is closer to the plane or below the plane
         now find the index of the edge up-lp in the edge list of lp */
627
      *ls = voronoi_get_edgeindex(c, (*up), (*us));
628

629
630
      /* if lp is also above the plane, replace up by lp and repeat the process
         until lp is below the plane */
631
      safewhile((*lw) == 1) {
632
        /* PATH 1.4 */
633
        *u = (*l);
634
635
636
637
638
639
640
641
642
        *up = (*lp);
        *us = 0;
        /* no while: PATH 1.4.0 */
        /* somewhere in while: PATH 1.4.1 */
        /* last valid option of while: PATH 1.4.2 */
        safewhile((*us) < (*ls) && (*l) >= (*u)) {
          *lp = voronoi_get_edge(c, (*up), (*us));
          *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l,
                                    teststack, &teststack_size);
643
          ++(*us);
644
645
        }
        if ((*l) >= (*u)) {
646
          ++(*us);
647
648
649
650
651
652
653
          /* no while: PATH 1.4.3 */
          /* somewhere in while: PATH 1.4.4 */
          /* last valid option of while: PATH 1.4.5 */
          safewhile((*us) < c->orders[(*up)] && (*l) >= (*u)) {
            *lp = voronoi_get_edge(c, (*up), (*us));
            *lw = voronoi_test_vertex(&c->vertices[3 * (*lp)], dx, r2, l,
                                      teststack, &teststack_size);
654
            ++(*us);
655
656
657
658
659
660
661
662
663
664
          }
          if ((*l) >= (*u)) {
            /* PATH 1.4.6 */
            message(
                "Cell completely gone! This should not happen. (l >= u, l = "
                "%g, u = %g)",
                (*l), (*u));
            return -1;
          }
        }
665
        --(*us);
666
667
        *ls = voronoi_get_edgeindex(c, (*up), (*us));
      }
668
669
      /* if lp is too close to the plane, replace up by lp and proceed to
         complicated setup */
670
671
672
673
674
675
      if ((*lw) == 0) {
        /* PATH 1.5 */
        *up = (*lp);
        complicated = 1;
      }
    } else { /* if(uw == 1) */
676

677
      /* PATH 2 */
678

679
680
681
      /* below: try to find a vertex above
         we test all edges of the current vertex stored in up (vertex 0) until
         we either find one above the plane or closer to the plane */
682

683
684
685
686
687
688
689
690
691
692
693
      *qp = voronoi_get_edge(c, (*up), 0);
      *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q, teststack,
                                &teststack_size);
      *us = 1;
      /* not in while: PATH 2.0 */
      /* somewhere in while: PATH 2.1 */
      /* last valid option of while: PATH 2.2 */
      safewhile((*us) < c->orders[(*up)] && (*u) >= (*q)) {
        *qp = voronoi_get_edge(c, (*up), (*us));
        *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q, teststack,
                                  &teststack_size);
694
        ++(*us);
695
696
697
      }
      if ((*u) >= (*q)) {
        /* PATH 2.3 */
698
699
700
701
702
        /* up is the closest vertex to the plane and is below the plane
           since the cell is convex, up is the closest vertex of all vertices of
           the cell
           this means that the entire cell is below the plane
           The cell is unaltered. */
703
704
        return 0;
      } else {
705
        /* the last increase in the loop pushed us too far, correct this */
706
        --(*us);
707
      }
708

709
      /* repeat the above process until qp is closer or above the plane */
710
711
712
713
714
715
716
717
718
719
720
721
722
      safewhile((*qw) == -1) {
        /* PATH 2.4 */
        *qs = voronoi_get_edgeindex(c, (*up), (*us));
        *u = (*q);
        *up = (*qp);
        *us = 0;
        /* no while: PATH 2.4.0 */
        /* somewhere in while: PATH 2.4.1 */
        /* last valid option of while: 2.4.2 */
        safewhile((*us) < (*qs) && (*u) >= (*q)) {
          *qp = voronoi_get_edge(c, (*up), (*us));
          *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q,
                                    teststack, &teststack_size);
723
          ++(*us);
724
725
        }
        if ((*u) >= (*q)) {
726
          ++(*us);
727
728
729
730
731
732
733
          /* no while: PATH 2.4.3 */
          /* somewhere in while: PATH 2.4.4 */
          /* last valid option of while: PATH 2.4.5 */
          safewhile((*us) < c->orders[(*up)] && (*u) >= (*q)) {
            *qp = voronoi_get_edge(c, (*up), (*us));
            *qw = voronoi_test_vertex(&c->vertices[3 * (*qp)], dx, r2, q,
                                      teststack, &teststack_size);
734
            ++(*us);
735
736
737
738
739
740
741
          }
          if ((*u) >= (*q)) {
            /* PATH 2.4.6 */
            /* cell unaltered */
            return 0;
          }
        }
742
        --(*us);
743
744
      }
      if ((*qw) == 1) {
745
        /* qp is above the plane: initialize lp to up and replace up by qp */
746
747
748
749
750
751
752
753
        *lp = (*up);
        *ls = (*us);
        *l = (*u);
        *up = (*qp);
        *us = voronoi_get_edgeindex(c, (*lp), (*ls));
        *u = (*q);
      } else {
        /* PATH 2.5 */
754
        /* too close to call: go to complicated setup */
755
756
757
        *up = (*qp);
        complicated = 1;
      }
758

759
    } /* if(uw == 1) */
760

761
  } /* if(uw == 0) */
762

763
764
765
766
767
  if (complicated) {
    return 2;
  } else {
    return 1;
  }
768
769
}

770
/**
771
772
 * @brief Intersect the given cell with the midplane between the cell generator
 * and a neighbouring cell at the given relative position and with the given ID.
773
 *
774
775
776
 * This method is the core of the Voronoi algorithm. If anything goes wrong
 * geometrically, it most likely goes wrong somewhere within this method.
 *
777
 * @param c 3D Voronoi cell.
778
779
780
781
 * @param odx The original relative distance vector between the cell generator
 * and the intersecting neighbour, as it is passed on to runner_iact_density
 * (remember: odx = pi->x - pj->x).
 * @param ngb ID of the intersecting neighbour (pj->id in runner_iact_density).
782
783
 */
__attribute__((always_inline)) INLINE void voronoi_intersect(
784
    struct voronoi_cell *c, const float *odx, unsigned long long ngb) {
785

786
787
  /* vector pointing from pi to the midpoint of the line segment between pi and
     pj. This corresponds to -0.5*odx */
788
  float dx[3];
789
  /* squared norm of dx */
790
  float r2;
791
792
793
794
  /* u: distance between the plane and the closest vertex above the plane (up)
     l: distance between the plane and the closest vertex below the plane (lp)
     q: distance between the plane and the vertex that is currently being
     tested (qp) */
795
  float u = 0.0f, l = 0.0f, q = 0.0f;
796
797
798
799
  /* up: index of the closest vertex above the plane
     us: index of the edge of vertex up that intersects the plane
     uw: result of the last orientation test involving vertex u
     same naming used for vertex l and vertex q */
800
801
  int up = -1, us = -1, uw = -1, lp = -1, ls = -1, lw = -1, qp = -1, qs = -1,
      qw = -1;
802
  /* auxiliary flag used to capture degeneracies */
803
804
  int complicated = -1;

805
806
  /* stack to store all vertices that have already been tested (debugging
     only) */
807
  float teststack[2 * VORONOI3D_MAXNUMVERT];
808
  /* size of the used part of the stack */
809
810
  int teststack_size = 0;

811
812
813
814
815
#ifdef VORONOI3D_EXPENSIVE_CHECKS
  voronoi_check_cell_consistency(c);
#endif

  /* initialize dx and r2 */
816
817
818
819
820
  dx[0] = -0.5f * odx[0];
  dx[1] = -0.5f * odx[1];
  dx[2] = -0.5f * odx[2];
  r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];

821
  /* find an intersected edge of the cell */
822
823
824
  int result = voronoi_intersect_find_closest_vertex(
      c, dx, r2, &u, &up, &us, &uw, &l, &lp, &ls, &lw, &q, &qp, &qs, &qw);
  if (result < 0) {
825
826
827
828
829
    /* the closest_vertex test only found vertices above the intersecting plane
       this would mean that the entire cell lies above the midplane of the line
       segment connecting a point inside the cell (the generator) and a point
       that could be inside or outside the cell (the neighbour). This is
       geometrically absurd and should NEVER happen. */
830
    voronoi_print_gnuplot_c(c);
831
832
    error("Error while searching intersected edge!");
  }
833
  if (result == 0) {
834
835
836
837
    /* no intersection */
    return;
  }
  if (result == 2) {
838
839
    complicated = 1;
  } else {
840
841
    complicated = 0;
  }
842

843
844
845
846
847
848
849
850
851
852
853
  /* At this point:
      up contains a vertex above the plane
      lp contains a vertex below the plane
      us and ls contain the index of the edge that connects up and lp, this edge
      is intersected by the midplane
      u and l contain the projected distances of up and lp to the midplane,
      along dx
     IF complicated is 1, up contains a vertex that is considered to be on the
     plane. All other variables can be considered to be uninitialized in this
     case. */

854
855
  int vindex = -1;
  int visitflags[VORONOI3D_MAXNUMVERT];
856
  int dstack[2 * VORONOI3D_MAXNUMVERT];
857
  int dstack_size = 0;
858
859
860
861
862
863
  float r = 0.0f;
  int cs = -1, rp = -1;
  int double_edge = 0;
  int i = -1, j = -1, k = -1;

  /* initialize visitflags */
864
  for (i = 0; i < VORONOI3D_MAXNUMVERT; ++i) {
865
866
867
868
869
    visitflags[i] = 0;
  }

  if (complicated) {

870
871
872
873
874
875
    /* We've entered the complicated setup, which means that somewhere along the
       way we found a vertex that is on or very close to the midplane. The index
       of that vertex is stored in up, all other variables are meaningless at
       this point. */

    /* first of all, we need to find a vertex which has edges that extend below
876
877
878
       the plane (since the remainder of our algorithm depends on that). This is
       not necessarily the case: in principle a vertex can only have edges that
       extend inside or above the plane.
879
880
881
882
       we create a stack of vertices to test (we use dstack for this), and add
       vertex up. For each vertex on the stack, we then traverse its edges. If
       the edge extends above the plane, we ignore it. If it extends below, we
       stop. If the edge lies in the plane, we add the vertex on the other end
883
       to the stack.
884
885
       We make sure that up contains the index of a vertex extending beyond the
       plane on exit. */
886
887
888
    dstack[dstack_size] = up;
    ++dstack_size;
    lw = 0;
889
    j = 0;
890
891
892
893
894
895
    safewhile(j < dstack_size && lw != -1) {
      up = dstack[j];
      for (i = 0; i < c->orders[up]; ++i) {
        lp = voronoi_get_edge(c, up, i);
        lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                                 &teststack_size);
896
897
898
        if (lw == -1) {
          /* jump out of the for loop */
          break;
899
900
        }
        if (lw == 0) {
901
          /* only add each vertex to the stack once */
902
903
904
905
906
907
908
          k = 0;
          safewhile(k < dstack_size && dstack[k] != lp) { ++k; }
          if (k == dstack_size) {
            dstack[dstack_size] = lp;
            ++dstack_size;
          }
        }
909
      }
910
      ++j;
911
    }
912
913

    /* we increased j after lw was calculated, so only the value of lw should be
914
       used to determine whether or not the loop was successful */
915
    if (lw != -1) {
916
917
918
919
920
921
922
923
924
      /* we did not find an edge that extends below the plane. There are two
         possible reasons for this: either all vertices of the cell lie above
         or inside the midplane of the segment connecting a point inside the
         cell (the generator) with a point inside or outside the cell (the
         neighbour). This is geometrically absurd.
         Another reason might be that somehow all vertices in the midplane only
         have edges that extend outwards. This is contradictory to the fact that
         a Voronoi cell is convex, and therefore also unacceptable.
         We conclude that we should NEVER end up here. */
925
      voronoi_print_cell(c);
926
      error("Unable to find a vertex below the midplane!");
927
    }
928
929
    /* reset the delete stack, we need it later on */
    dstack_size = 0;
930

931
    /* the search routine detected a vertex very close to or in the midplane
932
933
       the index of this vertex is stored in up
       we proceed by checking the edges of this vertex */
934

935
936
937
    lp = voronoi_get_edge(c, up, 0);
    lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                             &teststack_size);
938

939
    /* the first edge can be below, above or on the plane */
940
    if (lw != -1) {
941

942
      /* above or on the plane: we try to find one below the plane */
943

944
945
946
947
948
949
      rp = lw;
      i = 1;
      lp = voronoi_get_edge(c, up, i);
      lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                               &teststack_size);
      safewhile(lw != -1) {
950
        ++i;
951
        if (i == c->orders[up]) {
952
953
          /* none of the edges of up is below the plane. Since the cell is
             supposed to be convex, this means the entire cell is above or on
954
955
956
             the plane. This should not happen...
             Furthermore, we should really NEVER end up here, as in this case
             an error should already have be thrown above. */
957
958
959
960
961
962
963
964
965
966
967
968
969
          voronoi_print_gnuplot_c(c);
          error(
              "Cell completely gone! This should not happen. (i == "
              "c->order[up], i = %d, c->orders[up] = %d, up = %d)\n"
              "dx: [%g %g %g]\nv[up]: [%g %g %g]\nx: [%g %g %g]",
              i, c->orders[up], up, dx[0], dx[1], dx[2], c->vertices[3 * up],
              c->vertices[3 * up + 1], c->vertices[3 * up + 2], c->x[0],
              c->x[1], c->x[2]);
        }
        lp = voronoi_get_edge(c, up, i);
        lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                                 &teststack_size);
      }
970

971
972
973
974
975
976
977
978
      /* lp, l and lw now contain values corresponding to an edge below the
         plane
         rp contains the result of test_vertex for the first edge of up, for
         reference */

      /* we go on to the next edge of up, and see if we can find an edge that
         does not extend below the plane */

979
      j = i + 1;
980
      safewhile(j < c->orders[up] && lw == -1) {
981
982
983
        lp = voronoi_get_edge(c, up, j);
        lw = voronoi_test_vertex(&c->vertices[3 * lp], dx, r2, &l, teststack,
                                 &teststack_size);
984
        ++j;
985
      }
986

987
988
989
990
991
992
993
994
995
996
997
998
      if (lw != -1) {
        /* the last iteration increased j by 1 too many, correct this */
        --j;
      }

      /* j-i now contains the number of edges below the plane. We will replace
         up by a new vertex of order this number + 2 (since 2 new edges will be
         created inside the plane)
         however, we do not do this if there is exactly one edge that lies in
         the plane, and all other edges lie below, because in this case we can
         just keep vertex up as is */

999
      if (j == c->orders[up] && i == 1 && rp == 0) {
1000
        /* keep the order of up, and flag this event for later reference */
For faster browsing, not all history is shown. View entire blame