faunus
radical_tessellation_contact_construction.h
1 #ifndef VORONOTALT_RADICAL_TESSELLATION_CONTACT_CONSTRUCTION_H_
2 #define VORONOTALT_RADICAL_TESSELLATION_CONTACT_CONSTRUCTION_H_
3 
4 #include <vector>
5 
6 #include "basic_types_and_functions.h"
7 
8 namespace voronotalt
9 {
10 
11 class RadicalTessellationContactConstruction
12 {
13 public:
14  struct ContourPoint
15  {
16  SimplePoint p;
17  Float angle;
18  UnsignedInt left_id;
19  UnsignedInt right_id;
20  int indicator;
21 
22 
23  ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) noexcept :
24  p(p),
25  angle(FLOATCONST(0.0)),
26  left_id(left_id),
27  right_id(right_id),
28  indicator(0)
29  {
30  }
31  };
32 
33  typedef std::vector<ContourPoint> Contour;
34 
35  struct ContactDescriptor
36  {
37  Contour contour;
38  SimpleSphere intersection_circle_sphere;
39  SimplePoint intersection_circle_axis;
40  SimplePoint contour_barycenter;
41  Float sum_of_arc_angles;
42  Float area;
43  Float solid_angle_a;
44  Float solid_angle_b;
45  Float pyramid_volume_a;
46  Float pyramid_volume_b;
47  Float distance;
48  UnsignedInt flags;
49  UnsignedInt id_a;
50  UnsignedInt id_b;
51 
52  ContactDescriptor() noexcept :
53  sum_of_arc_angles(FLOATCONST(0.0)),
54  area(FLOATCONST(0.0)),
55  solid_angle_a(FLOATCONST(0.0)),
56  solid_angle_b(FLOATCONST(0.0)),
57  pyramid_volume_a(FLOATCONST(0.0)),
58  pyramid_volume_b(FLOATCONST(0.0)),
59  distance(FLOATCONST(0.0)),
60  flags(0),
61  id_a(0),
62  id_b(0)
63  {
64  }
65 
66  void clear() noexcept
67  {
68  id_a=0;
69  id_b=0;
70  contour.clear();
71  sum_of_arc_angles=FLOATCONST(0.0);
72  area=FLOATCONST(0.0);
73  solid_angle_a=FLOATCONST(0.0);
74  solid_angle_b=FLOATCONST(0.0);
75  pyramid_volume_a=FLOATCONST(0.0);
76  pyramid_volume_b=FLOATCONST(0.0);
77  distance=FLOATCONST(0.0);
78  flags=0;
79  }
80  };
81 
82  struct ContactDescriptorGraphics
83  {
84  std::vector<SimplePoint> outer_points;
85  SimplePoint barycenter;
86  SimplePoint plane_normal;
87 
88  ContactDescriptorGraphics() noexcept
89  {
90  }
91 
92  void clear() noexcept
93  {
94  outer_points.clear();
95  }
96  };
97 
98  static bool construct_contact_descriptor(
99  const std::vector<SimpleSphere>& spheres,
100  const std::vector<int>& spheres_exclusion_statuses,
101  const UnsignedInt a_id,
102  const UnsignedInt b_id,
103  const std::vector<ValuedID>& a_neighbor_collisions,
104  const Float max_circle_radius_restriction,
105  ContactDescriptor& result_contact_descriptor) noexcept
106  {
107  result_contact_descriptor.clear();
108  if(a_id<spheres.size() && b_id<spheres.size())
109  {
110  result_contact_descriptor.id_a=a_id;
111  result_contact_descriptor.id_b=b_id;
112  const SimpleSphere& a=spheres[a_id];
113  const SimpleSphere& b=spheres[b_id];
114  if(sphere_intersects_sphere(a, b) && !sphere_contains_sphere(a, b) && !sphere_contains_sphere(b, a))
115  {
116  result_contact_descriptor.intersection_circle_sphere=intersection_circle_of_two_spheres(a, b);
117  if(max_circle_radius_restriction>FLOATCONST(0.0))
118  {
119  result_contact_descriptor.intersection_circle_sphere.r=std::min(result_contact_descriptor.intersection_circle_sphere.r, max_circle_radius_restriction);
120  }
121  if(result_contact_descriptor.intersection_circle_sphere.r>FLOATCONST(0.0))
122  {
123  bool discarded=false;
124  bool contour_initialized=false;
125  {
126  bool nostop=true;
127  for(UnsignedInt i=0;i<a_neighbor_collisions.size() && !discarded && nostop;i++)
128  {
129  const UnsignedInt neighbor_id=a_neighbor_collisions[i].index;
130  if(neighbor_id!=b_id && (neighbor_id>=spheres_exclusion_statuses.size() || spheres_exclusion_statuses[neighbor_id]==0))
131  {
132  const SimpleSphere& c=spheres[neighbor_id];
133  if(sphere_intersects_sphere(result_contact_descriptor.intersection_circle_sphere, c) && sphere_intersects_sphere(b, c))
134  {
135  if(sphere_contains_sphere(c, a) || sphere_contains_sphere(c, b))
136  {
137  discarded=true;
138  }
139  else
140  {
141  const SimplePoint neighbor_ac_plane_center=center_of_intersection_circle_of_two_spheres(a, c);
142  const SimplePoint neighbor_ac_plane_normal=unit_point(sub_of_points(c.p, a.p));
143  const Float cos_val=dot_product(unit_point(sub_of_points(result_contact_descriptor.intersection_circle_sphere.p, a.p)), unit_point(sub_of_points(neighbor_ac_plane_center, a.p)));
144  if(std::abs(cos_val)<FLOATCONST(1.0))
145  {
146  const Float l=std::abs(signed_distance_from_point_to_plane(neighbor_ac_plane_center, neighbor_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p));
147  const Float xl=l/std::sqrt(1-(cos_val*cos_val));
148  if(xl>=result_contact_descriptor.intersection_circle_sphere.r)
149  {
150  if(halfspace_of_point(neighbor_ac_plane_center, neighbor_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p)>=0)
151  {
152  discarded=true;
153  }
154  }
155  else
156  {
157  if(!contour_initialized)
158  {
159  result_contact_descriptor.intersection_circle_axis=unit_point(sub_of_points(b.p, a.p));
160  init_contour_from_base_and_axis(a_id, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.intersection_circle_axis, result_contact_descriptor.contour);
161  contour_initialized=true;
162  }
163  else
164  {
165  nostop=(test_if_contour_is_still_cuttable(a.p, neighbor_ac_plane_center, result_contact_descriptor.contour));
166  }
167  if(nostop)
168  {
169  mark_and_cut_contour(neighbor_ac_plane_center, neighbor_ac_plane_normal, neighbor_id, result_contact_descriptor.contour);
170  if(contour_initialized && result_contact_descriptor.contour.empty())
171  {
172  discarded=true;
173  }
174  }
175  }
176  }
177  else
178  {
179  if(halfspace_of_point(neighbor_ac_plane_center, neighbor_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p)>0)
180  {
181  discarded=true;
182  }
183  }
184  }
185  }
186  }
187  }
188  }
189  if(!discarded)
190  {
191  if(!contour_initialized)
192  {
193  result_contact_descriptor.intersection_circle_axis=unit_point(sub_of_points(b.p, a.p));
194  result_contact_descriptor.contour_barycenter=result_contact_descriptor.intersection_circle_sphere.p;
195  result_contact_descriptor.sum_of_arc_angles=(PIVALUE*FLOATCONST(2.0));
196  result_contact_descriptor.area=result_contact_descriptor.intersection_circle_sphere.r*result_contact_descriptor.intersection_circle_sphere.r*PIVALUE;
197  }
198  else
199  {
200  if(!result_contact_descriptor.contour.empty())
201  {
202  restrict_contour_to_circle(result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.intersection_circle_axis, a_id, result_contact_descriptor.contour, result_contact_descriptor.sum_of_arc_angles);
203  if(!result_contact_descriptor.contour.empty())
204  {
205  result_contact_descriptor.area=calculate_contour_area(result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour, result_contact_descriptor.contour_barycenter);
206  }
207  }
208  }
209 
210  if(result_contact_descriptor.area>FLOATCONST(0.0))
211  {
212  result_contact_descriptor.solid_angle_a=calculate_contour_solid_angle(a, b, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour);
213  result_contact_descriptor.solid_angle_b=calculate_contour_solid_angle(b, a, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour);
214  result_contact_descriptor.pyramid_volume_a=distance_from_point_to_point(a.p, result_contact_descriptor.intersection_circle_sphere.p)*result_contact_descriptor.area/FLOATCONST(3.0)*(result_contact_descriptor.solid_angle_a<FLOATCONST(0.0) ? FLOATCONST(-1.0) : FLOATCONST(1.0));
215  result_contact_descriptor.pyramid_volume_b=distance_from_point_to_point(b.p, result_contact_descriptor.intersection_circle_sphere.p)*result_contact_descriptor.area/FLOATCONST(3.0)*(result_contact_descriptor.solid_angle_b<FLOATCONST(0.0) ? FLOATCONST(-1.0) : FLOATCONST(1.0));
216  result_contact_descriptor.flags=(check_if_contour_is_central(result_contact_descriptor.intersection_circle_sphere.p, result_contact_descriptor.contour, result_contact_descriptor.contour_barycenter) ? 1 : 0);
217  }
218 
219  result_contact_descriptor.distance=distance_from_point_to_point(a.p, b.p);
220  }
221  }
222  }
223  }
224  return (result_contact_descriptor.area>FLOATCONST(0.0));
225  }
226 
227  static bool construct_contact_descriptor_graphics(const ContactDescriptor& contact_descriptor, const Float length_step, ContactDescriptorGraphics& result_contact_descriptor_graphics) noexcept
228  {
229  result_contact_descriptor_graphics.clear();
230  if(contact_descriptor.area>FLOATCONST(0.0))
231  {
232  result_contact_descriptor_graphics.plane_normal=contact_descriptor.intersection_circle_axis;
233  const Float angle_step=std::max(std::min(length_step/contact_descriptor.intersection_circle_sphere.r, PIVALUE/FLOATCONST(3.0)), PIVALUE/FLOATCONST(36.0));
234  if(contact_descriptor.contour.empty())
235  {
236  const SimplePoint first_point=point_and_number_product(any_normal_of_vector(contact_descriptor.intersection_circle_axis), contact_descriptor.intersection_circle_sphere.r);
237  result_contact_descriptor_graphics.outer_points.reserve(static_cast<int>((PIVALUE*FLOATCONST(2.0))/angle_step)+2);
238  result_contact_descriptor_graphics.outer_points.push_back(sum_of_points(contact_descriptor.intersection_circle_sphere.p, first_point));
239  for(Float rotation_angle=angle_step;rotation_angle<(PIVALUE*FLOATCONST(2.0));rotation_angle+=angle_step)
240  {
241  result_contact_descriptor_graphics.outer_points.push_back(sum_of_points(contact_descriptor.intersection_circle_sphere.p, rotate_point_around_axis(contact_descriptor.intersection_circle_axis, rotation_angle, first_point)));
242  }
243  result_contact_descriptor_graphics.barycenter=contact_descriptor.intersection_circle_sphere.p;
244  }
245  else
246  {
247  if(contact_descriptor.sum_of_arc_angles>FLOATCONST(0.0))
248  {
249  result_contact_descriptor_graphics.outer_points.reserve(static_cast<UnsignedInt>(contact_descriptor.sum_of_arc_angles/angle_step)+contact_descriptor.contour.size()+4);
250  }
251  else
252  {
253  result_contact_descriptor_graphics.outer_points.reserve(contact_descriptor.contour.size());
254  }
255  for(UnsignedInt i=0;i<contact_descriptor.contour.size();i++)
256  {
257  const ContourPoint& pr=contact_descriptor.contour[i];
258  result_contact_descriptor_graphics.outer_points.push_back(pr.p);
259  if(pr.angle>FLOATCONST(0.0))
260  {
261  if(pr.angle>angle_step)
262  {
263  const SimplePoint first_v=sub_of_points(pr.p, contact_descriptor.intersection_circle_sphere.p);
264  for(Float rotation_angle=angle_step;rotation_angle<pr.angle;rotation_angle+=angle_step)
265  {
266  result_contact_descriptor_graphics.outer_points.push_back(sum_of_points(contact_descriptor.intersection_circle_sphere.p, rotate_point_around_axis(contact_descriptor.intersection_circle_axis, rotation_angle, first_v)));
267  }
268  }
269  }
270  }
271  result_contact_descriptor_graphics.barycenter.x=FLOATCONST(0.0);
272  result_contact_descriptor_graphics.barycenter.y=FLOATCONST(0.0);
273  result_contact_descriptor_graphics.barycenter.z=FLOATCONST(0.0);
274  if(!result_contact_descriptor_graphics.outer_points.empty())
275  {
276  for(UnsignedInt i=0;i<result_contact_descriptor_graphics.outer_points.size();i++)
277  {
278  set_close_to_equal(result_contact_descriptor_graphics.outer_points[i].x, FLOATCONST(0.0));
279  set_close_to_equal(result_contact_descriptor_graphics.outer_points[i].y, FLOATCONST(0.0));
280  set_close_to_equal(result_contact_descriptor_graphics.outer_points[i].z, FLOATCONST(0.0));
281  result_contact_descriptor_graphics.barycenter.x+=result_contact_descriptor_graphics.outer_points[i].x;
282  result_contact_descriptor_graphics.barycenter.y+=result_contact_descriptor_graphics.outer_points[i].y;
283  result_contact_descriptor_graphics.barycenter.z+=result_contact_descriptor_graphics.outer_points[i].z;
284  }
285  result_contact_descriptor_graphics.barycenter.x/=static_cast<Float>(result_contact_descriptor_graphics.outer_points.size());
286  result_contact_descriptor_graphics.barycenter.y/=static_cast<Float>(result_contact_descriptor_graphics.outer_points.size());
287  result_contact_descriptor_graphics.barycenter.z/=static_cast<Float>(result_contact_descriptor_graphics.outer_points.size());
288  }
289  }
290  }
291  return (!result_contact_descriptor_graphics.outer_points.empty());
292  }
293 
294 private:
295  static void init_contour_from_base_and_axis(
296  const UnsignedInt a_id,
297  const SimpleSphere& base,
298  const SimplePoint& axis,
299  Contour& result) noexcept
300  {
301  const SimplePoint first_point=point_and_number_product(any_normal_of_vector(axis), base.r*FLOATCONST(1.19));
302  const Float angle_step=PIVALUE/FLOATCONST(3.0);
303  result.push_back(ContourPoint(sum_of_points(base.p, first_point), a_id, a_id));
304  for(Float rotation_angle=angle_step;rotation_angle<(PIVALUE*FLOATCONST(2.0));rotation_angle+=angle_step)
305  {
306  result.push_back(ContourPoint(sum_of_points(base.p, rotate_point_around_axis(axis, rotation_angle, first_point)), a_id, a_id));
307  }
308  }
309 
310  static bool mark_and_cut_contour(
311  const SimplePoint& ac_plane_center,
312  const SimplePoint& ac_plane_normal,
313  const UnsignedInt c_id,
314  Contour& contour) noexcept
315  {
316  const UnsignedInt outsiders_count=mark_contour(ac_plane_center, ac_plane_normal, c_id, contour);
317  if(outsiders_count>0)
318  {
319  if(outsiders_count<contour.size())
320  {
321  cut_contour(ac_plane_center, ac_plane_normal, c_id, contour);
322  }
323  else
324  {
325  contour.clear();
326  }
327  return true;
328  }
329  return false;
330  }
331 
332  static UnsignedInt mark_contour(
333  const SimplePoint& ac_plane_center,
334  const SimplePoint& ac_plane_normal,
335  const UnsignedInt c_id,
336  Contour& contour) noexcept
337  {
338  UnsignedInt outsiders_count=0;
339  for(Contour::iterator it=contour.begin();it!=contour.end();++it)
340  {
341  if(halfspace_of_point(ac_plane_center, ac_plane_normal, it->p)>=0)
342  {
343  it->left_id=c_id;
344  it->right_id=c_id;
345  outsiders_count++;
346  }
347  }
348  return outsiders_count;
349  }
350 
351  static void cut_contour(
352  const SimplePoint& ac_plane_center,
353  const SimplePoint& ac_plane_normal,
354  const UnsignedInt c_id,
355  Contour& contour) noexcept
356  {
357  if(contour.size()<3)
358  {
359  return;
360  }
361 
362  UnsignedInt i_start=0;
363  while(i_start<contour.size() && !(contour[i_start].left_id==c_id && contour[i_start].right_id==c_id))
364  {
365  i_start++;
366  }
367 
368  if(i_start>=contour.size())
369  {
370  return;
371  }
372 
373  UnsignedInt i_end=(contour.size()-1);
374  while(i_end>0 && !(contour[i_end].left_id==c_id && contour[i_end].right_id==c_id))
375  {
376  i_end--;
377  }
378 
379  if(i_start==0 && i_end==(contour.size()-1))
380  {
381  i_end=0;
382  while((i_end+1)<contour.size() && contour[i_end+1].left_id==c_id && contour[i_end+1].right_id==c_id)
383  {
384  i_end++;
385  }
386 
387  i_start=(contour.size()-1);
388  while(i_start>0 && contour[i_start-1].left_id==c_id && contour[i_start-1].right_id==c_id)
389  {
390  i_start--;
391  }
392  }
393 
394  if(i_start==i_end)
395  {
396  contour.insert(contour.begin()+i_start, contour[i_start]);
397  i_end=i_start+1;
398  }
399  else if(i_start<i_end)
400  {
401  if(i_start+1<i_end)
402  {
403  contour.erase(contour.begin()+i_start+1, contour.begin()+i_end);
404  }
405  i_end=i_start+1;
406  }
407  else if(i_start>i_end)
408  {
409  if(i_start+1<contour.size())
410  {
411  contour.erase(contour.begin()+i_start+1, contour.end());
412  }
413  if(i_end>0)
414  {
415  contour.erase(contour.begin(), contour.begin()+i_end);
416 
417  }
418  i_start=contour.size()-1;
419  i_end=0;
420  }
421 
422  {
423  const UnsignedInt i_left=((i_start)>0 ? (i_start-1) : (contour.size()-1));
424  contour[i_start]=ContourPoint(intersection_of_plane_and_segment(ac_plane_center, ac_plane_normal, contour[i_start].p, contour[i_left].p), contour[i_left].right_id, contour[i_start].left_id);
425  }
426 
427  {
428  const UnsignedInt i_right=((i_end+1)<contour.size() ? (i_end+1) : 0);
429  contour[i_end]=ContourPoint(intersection_of_plane_and_segment(ac_plane_center, ac_plane_normal, contour[i_end].p, contour[i_right].p), contour[i_end].right_id, contour[i_right].left_id);
430  }
431 
432  if(!greater(squared_distance_from_point_to_point(contour[i_start].p, contour[i_end].p), FLOATCONST(0.0)))
433  {
434  contour[i_start].right_id=contour[i_end].right_id;
435  contour.erase(contour.begin()+i_end);
436  }
437  }
438 
439  static bool test_if_contour_is_still_cuttable(const SimplePoint& a_center, const SimplePoint& closest_possible_cut_point, const Contour& contour) noexcept
440  {
441  bool cuttable=false;
442  const Float dist_threshold=squared_distance_from_point_to_point(a_center, closest_possible_cut_point);
443  for(UnsignedInt i=0;!cuttable && i<contour.size();i++)
444  {
445  cuttable=squared_distance_from_point_to_point(a_center, contour[i].p)>=dist_threshold;
446  }
447  return cuttable;
448  }
449 
450  static bool restrict_contour_to_circle(
451  const SimpleSphere& ic_sphere,
452  const SimplePoint& ic_axis,
453  const UnsignedInt a_id,
454  Contour& contour,
455  Float& sum_of_arc_angles) noexcept
456  {
457  UnsignedInt outsiders_count=0;
458  for(UnsignedInt i=0;i<contour.size();i++)
459  {
460  if(squared_distance_from_point_to_point(contour[i].p, ic_sphere.p)<=(ic_sphere.r*ic_sphere.r))
461  {
462  contour[i].indicator=0;
463  }
464  else
465  {
466  contour[i].indicator=1;
467  outsiders_count++;
468  }
469  }
470 
471  if(outsiders_count>0)
472  {
473  UnsignedInt insertions_count=0;
474  {
475  UnsignedInt i=0;
476  while(i<contour.size())
477  {
478  ContourPoint& pr1=contour[i];
479  ContourPoint& pr2=contour[(i+1)<contour.size() ? (i+1) : 0];
480  if(pr1.indicator==1 || pr2.indicator==1)
481  {
482  if(pr1.indicator==1 && pr2.indicator==1)
483  {
484  SimplePoint mp;
485  if(project_point_inside_line(ic_sphere.p, pr1.p, pr2.p, mp))
486  {
487  if(squared_distance_from_point_to_point(mp, ic_sphere.p)<=(ic_sphere.r*ic_sphere.r))
488  {
489  SimplePoint ip1;
490  SimplePoint ip2;
491  if(intersect_segment_with_circle(ic_sphere, mp, pr1.p, ip1) && intersect_segment_with_circle(ic_sphere, mp, pr2.p, ip2))
492  {
493  const UnsignedInt pr2_left_id=pr2.left_id;
494  contour.insert(((i+1)<contour.size() ? (contour.begin()+(i+1)) : contour.end()), 2, ContourPoint(ip1, a_id, pr1.right_id));
495  contour[i+2]=ContourPoint(ip2, pr2_left_id, a_id);
496  insertions_count+=2;
497  i+=2;
498  }
499  }
500  }
501  }
502  else
503  {
504  if(pr1.indicator==1 && pr2.indicator!=1)
505  {
506  SimplePoint ip;
507  if(intersect_segment_with_circle(ic_sphere, pr2.p, pr1.p, ip))
508  {
509  contour.insert(((i+1)<contour.size() ? (contour.begin()+(i+1)) : contour.end()), ContourPoint(ip, a_id, pr1.right_id));
510  insertions_count++;
511  i++;
512  }
513  else
514  {
515  pr2.left_id=a_id;
516  pr2.right_id=pr1.right_id;
517  }
518  }
519  else if(pr1.indicator!=1 && pr2.indicator==1)
520  {
521  SimplePoint ip;
522  if(intersect_segment_with_circle(ic_sphere, pr1.p, pr2.p, ip))
523  {
524  contour.insert(((i+1)<contour.size() ? (contour.begin()+(i+1)) : contour.end()), ContourPoint(ip, pr2.left_id, a_id));
525  insertions_count++;
526  i++;
527  }
528  else
529  {
530  pr1.left_id=pr2.left_id;
531  pr1.right_id=a_id;
532  }
533  }
534  }
535  }
536  i++;
537  }
538  }
539  if(insertions_count==0)
540  {
541  contour.clear();
542  }
543  else if(!contour.empty())
544  {
545  {
546  UnsignedInt i=(contour.size()-1);
547  while(i<contour.size())
548  {
549  if(contour[i].indicator==1)
550  {
551  contour.erase(contour.begin()+i);
552  }
553  i=(i>0 ? (i-1) : contour.size());
554  }
555  }
556  if(contour.size()<2)
557  {
558  contour.clear();
559  }
560  {
561  UnsignedInt i=0;
562  while(i<contour.size())
563  {
564  ContourPoint& pr1=contour[i];
565  ContourPoint& pr2=contour[(i+1)<contour.size() ? (i+1) : 0];
566  if(pr1.right_id==a_id && pr2.left_id==a_id)
567  {
568  pr1.angle=directed_angle(ic_sphere.p, pr1.p, pr2.p, sum_of_points(ic_sphere.p, ic_axis));
569  sum_of_arc_angles+=pr1.angle;
570  pr1.indicator=2;
571  pr2.indicator=3;
572  }
573  i++;
574  }
575  }
576  if(greater_or_equal(sum_of_arc_angles, PIVALUE*FLOATCONST(2.0)) || (contour.size()>2 && equal(sum_of_arc_angles, PIVALUE*FLOATCONST(2.0), 0.001)))
577  {
578  sum_of_arc_angles=FLOATCONST(2.0);
579  contour.clear();
580  }
581  }
582  }
583 
584  return (outsiders_count>0);
585  }
586 
587  static Float calculate_contour_area(const SimpleSphere& ic_sphere, const Contour& contour, SimplePoint& contour_barycenter) noexcept
588  {
589  Float area=FLOATCONST(0.0);
590 
591  contour_barycenter.x=FLOATCONST(0.0);
592  contour_barycenter.y=FLOATCONST(0.0);
593  contour_barycenter.z=FLOATCONST(0.0);
594  for(UnsignedInt i=0;i<contour.size();i++)
595  {
596  contour_barycenter.x+=contour[i].p.x;
597  contour_barycenter.y+=contour[i].p.y;
598  contour_barycenter.z+=contour[i].p.z;
599  }
600  contour_barycenter.x/=static_cast<Float>(contour.size());
601  contour_barycenter.y/=static_cast<Float>(contour.size());
602  contour_barycenter.z/=static_cast<Float>(contour.size());
603 
604  for(UnsignedInt i=0;i<contour.size();i++)
605  {
606  const ContourPoint& pr1=contour[i];
607  const ContourPoint& pr2=contour[(i+1)<contour.size() ? (i+1) : 0];
608  area+=triangle_area(contour_barycenter, pr1.p, pr2.p);
609  if(pr1.angle>FLOATCONST(0.0))
610  {
611  area+=ic_sphere.r*ic_sphere.r*(pr1.angle-std::sin(pr1.angle))*0.5;
612  }
613  }
614 
615  return area;
616  }
617 
618  static Float calculate_contour_solid_angle(const SimpleSphere& a, const SimpleSphere& b, const SimpleSphere& ic_sphere, const Contour& contour) noexcept
619  {
620  Float turn_angle=FLOATCONST(0.0);
621 
622  if(!contour.empty())
623  {
624  for(UnsignedInt i=0;i<contour.size();i++)
625  {
626  const ContourPoint& pr0=contour[(i>0) ? (i-1) : (contour.size()-1)];
627  const ContourPoint& pr1=contour[i];
628  const ContourPoint& pr2=contour[(i+1)<contour.size() ? (i+1) : 0];
629 
630  if(pr0.angle>FLOATCONST(0.0))
631  {
632  SimplePoint d=cross_product(sub_of_points(b.p, a.p), sub_of_points(pr1.p, ic_sphere.p));
633  if((pr0.angle<PIVALUE && dot_product(d, sub_of_points(pr0.p, pr1.p))<FLOATCONST(0.0)) || (pr0.angle>PIVALUE && dot_product(d, sub_of_points(pr0.p, pr1.p))>FLOATCONST(0.0)))
634  {
635  d=point_and_number_product(d, FLOATCONST(-1.0));
636  }
637  turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, sum_of_points(pr1.p, d), pr2.p));
638  }
639  else if(pr1.angle>FLOATCONST(0.0))
640  {
641  SimplePoint d=cross_product(sub_of_points(b.p, a.p), sub_of_points(pr1.p, ic_sphere.p));
642  if((pr1.angle<PIVALUE && dot_product(d, sub_of_points(pr2.p, pr1.p))<FLOATCONST(0.0)) || (pr1.angle>PIVALUE && dot_product(d, sub_of_points(pr2.p, pr1.p))>FLOATCONST(0.0)))
643  {
644  d=point_and_number_product(d, FLOATCONST(-1.0));
645  }
646  turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, pr0.p, sum_of_points(pr1.p, d)));
647  turn_angle+=pr1.angle*(distance_from_point_to_point(a.p, ic_sphere.p)/a.r);
648  }
649  else
650  {
651  turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, pr0.p, pr2.p));
652  }
653  }
654  }
655  else
656  {
657  turn_angle=(PIVALUE*FLOATCONST(2.0))*(distance_from_point_to_point(a.p, ic_sphere.p)/a.r);
658  }
659 
660  Float solid_angle=((PIVALUE*FLOATCONST(2.0))-turn_angle);
661 
662  if(dot_product(sub_of_points(ic_sphere.p, a.p), sub_of_points(ic_sphere.p, b.p))>FLOATCONST(0.0) && squared_distance_from_point_to_point(ic_sphere.p, a.p)<squared_distance_from_point_to_point(ic_sphere.p, b.p))
663  {
664  solid_angle=(FLOATCONST(0.0)-solid_angle);
665  }
666 
667  return solid_angle;
668  }
669 
670  static bool check_if_contour_is_central(const SimplePoint& center, const Contour& contour, const SimplePoint& contour_barycenter) noexcept
671  {
672  bool central=false;
673  for(UnsignedInt i=0;i<contour.size() && !central;i++)
674  {
675  central=central || greater(contour[i].angle, PIVALUE);
676  }
677  if(!central)
678  {
679  central=true;
680  for(UnsignedInt i=0;i<contour.size() && central;i++)
681  {
682  const UnsignedInt j=((i+1)<contour.size() ? (i+1) : 0);
683  const SimplePoint u_ij=unit_point(sub_of_points(contour[j].p, contour[i].p));
684  const SimplePoint n_ijb=sub_of_points(contour_barycenter, sum_of_points(contour[i].p, point_and_number_product(u_ij, dot_product(u_ij, sub_of_points(contour_barycenter, contour[i].p)))));
685  if(dot_product(n_ijb, sub_of_points(center, contour[i].p))<FLOATCONST(0.0))
686  {
687  central=false;
688  }
689  }
690  }
691  return central;
692  }
693 };
694 
695 }
696 
697 #endif /* VORONOTALT_RADICAL_TESSELLATION_CONTACT_CONSTRUCTION_H_ */
voronotalt::RadicalTessellationContactConstruction::ContactDescriptor
Definition: radical_tessellation_contact_construction.h:35
voronotalt::SimplePoint
Definition: basic_types_and_functions.h:34
voronotalt::RadicalTessellationContactConstruction::ContourPoint
Definition: radical_tessellation_contact_construction.h:14
voronotalt::RadicalTessellationContactConstruction
Definition: radical_tessellation_contact_construction.h:11
voronotalt::RadicalTessellationContactConstruction::ContactDescriptorGraphics
Definition: radical_tessellation_contact_construction.h:82
voronotalt
Definition: basic_types_and_functions.h:19
voronotalt::SimpleSphere
Definition: basic_types_and_functions.h:49
Generated by   doxygen 1.8.13