diff --git a/docs/_docs/analysis.md b/docs/_docs/analysis.md index 441224776..362032238 100644 --- a/docs/_docs/analysis.md +++ b/docs/_docs/analysis.md @@ -364,6 +364,22 @@ is filed to disk, `sasa_histogram.dat`. `file` | Optionally stream area for each `nstep` to file (`.dat|.dat.gz`) `radius=1.4` | Probe radius (Å) +### Voronoi Tessellation (experimental) + +Performs Voronoi tessellation to detect contacts and surface areas using +the [Voronota-LT library](https://doi.org/10/mq8k). +Only the total surface area is currently reported, but more options are planned for future +releases. +The algorithm is significantly faster than the above `sasa` analysis. + +`voronoi` | Description +-------------- | --------------------------- +`nstep` | Interval between samples +`nskip=0` | Number of initial steps excluded from the analysis +`file` | Optionally stream surface area for each `nstep` to file (`.dat|.dat.gz`) +`radius=1.4` | Probe radius (Å) + + ## Charge Properties ### Molecular Multipoles diff --git a/docs/schema.yml b/docs/schema.yml index bd0bedc55..ed67c83c9 100644 --- a/docs/schema.yml +++ b/docs/schema.yml @@ -1263,6 +1263,23 @@ properties: nskip: {type: integer, default: 0, description: Number of steps to initially skip} required: [dL, molecule, nstep] additionalProperties: false + + voronoi: + description: "Voronoi tessellation using Voronota" + type: object + properties: + nstep: {type: integer, description: "Interval between samples"} + nskip: {type: integer, default: 0, description: "Initial steps to skip"} + radius: + type: number + default: 1.4 + description: "Probe radius (Å)" + file: + type: string + pattern: "(.*?)\\.(dat|dat.gz)$" + description: "Streaming filename (.dat|.dat.gz)" + required: [nstep] + additionalProperties: false widom: description: "Widom ghost particle insertion" diff --git a/examples/sasa/sasa.out.json b/examples/sasa/sasa.out.json index a78cbdd2f..9aa5812d9 100644 --- a/examples/sasa/sasa.out.json +++ b/examples/sasa/sasa.out.json @@ -6,12 +6,21 @@ "nstep": 1, "radius": 2.0, "reference": "doi:10/dbjh", - "relative time": 0.592, "samples": 1, "slices_per_atom": 20, "⟨SASA²⟩-⟨SASA⟩²": 0.0, "⟨SASA⟩": 2177.145143881812 } + }, + { + "voronota": { + "nstep": 1, + "radius": 2.0, + "reference": "doi:10/mq8k", + "samples": 1, + "⟨SASA²⟩-⟨SASA⟩²": 0.0, + "⟨SASA⟩": 2186.5484868984963 + } } ], "compiler": "Apple LLVM 15.0.0 (clang-1500.0.40.1)", diff --git a/examples/sasa/sasa.yml b/examples/sasa/sasa.yml index 04c4c53e6..58f634b9c 100755 --- a/examples/sasa/sasa.yml +++ b/examples/sasa/sasa.yml @@ -26,4 +26,8 @@ analysis: molecule: mymolecule radius: 2.0 nstep: 1 + - voronoi: + file: voronoi_sasa.dat + radius: 2.0 + nstep: 1 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 62d244ed3..1a98c37bb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,6 +8,8 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) target_link_libraries(project_options INTERFACE ZLIB::ZLIB) target_link_libraries(project_options INTERFACE ${CMAKE_THREAD_LIBS_INIT}) +include_directories(SYSTEM ${CMAKE_SOURCE_DIR}/voronotalt) + # ========== pre-compiler system headers ========== option(ENABLE_PCH "Enable Precompiled Headers" OFF) @@ -60,7 +62,7 @@ set(objs actions.cpp analysis.cpp average.cpp atomdata.cpp auxiliary.cpp bonds.c chainmove.cpp clustermove.cpp core.cpp forcemove.cpp units.cpp energy.cpp externalpotential.cpp geometry.cpp group.cpp io.cpp molecule.cpp montecarlo.cpp move.cpp mpicontroller.cpp particle.cpp penalty.cpp potentials.cpp random.cpp reactioncoordinate.cpp regions.cpp rotate.cpp sasa.cpp - scatter.cpp smart_montecarlo.cpp space.cpp speciation.cpp spherocylinder.cpp tensor.cpp) + scatter.cpp smart_montecarlo.cpp space.cpp speciation.cpp spherocylinder.cpp tensor.cpp voronota.cpp) set(hdrs actions.h analysis.h average.h atomdata.h auxiliary.h bonds.h celllist.h celllistimpl.h chainmove.h clustermove.h core.h forcemove.h energy.h externalpotential.h geometry.h group.h io.h diff --git a/src/analysis.cpp b/src/analysis.cpp index 8269e74fa..d4e7f73b2 100644 --- a/src/analysis.cpp +++ b/src/analysis.cpp @@ -25,7 +25,6 @@ #include #include #include - #include #include #include @@ -194,6 +193,8 @@ std::unique_ptr createAnalysis(const std::string& name, const json& j, return std::make_unique(j, spc, pot); } else if (name == "virtualtranslate") { return std::make_unique(j, spc, pot); + } else if (name == "voronoi") { + return std::make_unique(j, spc); } else if (name == "widom") { return std::make_unique(j, spc, pot); } else if (name == "xtcfile") { @@ -1879,6 +1880,7 @@ SASAAnalysis::SASAAnalysis(const double probe_radius, const int slices_per_atom, break; } setPolicy(selected_policy); + cite = "doi:10/dbjh"; } SASAAnalysis::SASAAnalysis(const json& j, const Space& spc) @@ -1886,7 +1888,6 @@ SASAAnalysis::SASAAnalysis(const json& j, const Space& spc) j.value("policy", Policies::INVALID), spc) { from_json(j); policy->from_json(j); - cite = "doi:10/dbjh"; } void SASAAnalysis::_to_json(json& json_ouput) const { diff --git a/src/analysis.h b/src/analysis.h index 3465f2a3e..bd5dc0e25 100644 --- a/src/analysis.h +++ b/src/analysis.h @@ -262,7 +262,7 @@ class ElectricPotential : public Analysis { unsigned int calculations_per_sample_event = 1; std::string file_prefix; //!< Output filename prefix for potential histogram and correlation struct Target { - Point position; //!< Target position + Point position; //!< Target position Average mean_potential; //!< mean potential at position std::unique_ptr> potential_histogram; //!< Histogram of observed potentials }; @@ -1086,6 +1086,36 @@ class SavePenaltyEnergy : public Analysis { SavePenaltyEnergy(const json& j, const Space& spc, const Energy::Hamiltonian& pot); }; +/** + * @brief Analysis of Vorononoi tessellation using the Voronota-LT library + * + * @todo Currently only SASA information is reported in the output. Add more information! + * + * https://doi.org/10/mq8k + */ +class Voronota : public Analysis { + private: + struct Averages { + Average area; + Average area_squared; + }; //!< Placeholder class for average properties + Averages average_data; //!< Stores all averages for the selected molecule + + std::unique_ptr output_stream; //!< output stream + double probe_radius; //!< radius of the probe sphere + std::string filename; //!< output file name + bool use_pbc = false; //!< Is the cell periodic? + + void _to_json(json& j) const override; + void _from_json(const json& input) override; + void _to_disk() override; + void _sample() override; + + public: + Voronota(const json& j, const Space& spc); + Voronota(double probe_radius, const Space& spc); +}; + } // namespace analysis } // namespace Faunus diff --git a/src/voronota.cpp b/src/voronota.cpp new file mode 100644 index 000000000..515b09c6d --- /dev/null +++ b/src/voronota.cpp @@ -0,0 +1,85 @@ +#include "analysis.h" +#include "mpicontroller.h" +#include +#include + +namespace Faunus::analysis { + +Voronota::Voronota(double probe_radius, const Faunus::Space& spc) + : Analysis(spc, "voronota") + , probe_radius(probe_radius) { + cite = "doi:10/mq8k"; + auto n_pbc = spc.geometry.asSimpleGeometry()->boundary_conditions.isPeriodic().count(); + switch (n_pbc) { + case 0: + faunus_logger->debug("{}: No PBC detected", name); + use_pbc = false; + break; + case 3: + faunus_logger->debug("{}: 3D PBC detected", name); + use_pbc = true; + break; + default: + faunus_logger->warn("{}: Non-uniform PBC is currently ignored - be careful!", name); + use_pbc = false; + } +} + +Voronota::Voronota(const Faunus::json& input, const Faunus::Space& spc) + : Voronota(input.value("radius", 1.4_angstrom), spc) { + from_json(input); +} + +void Voronota::_from_json(const Faunus::json& input) { + if (filename = input.value("file", ""s); !filename.empty()) { + output_stream = IO::openCompressedOutputStream(MPI::prefix + filename); + *output_stream << "# step SASA\n"; + } +} + +void Voronota::_to_json(json& json_ouput) const { + if (!average_data.area.empty()) { + json_ouput = {{"⟨SASA⟩", average_data.area.avg()}, + {"⟨SASA²⟩-⟨SASA⟩²", average_data.area_squared.avg() - std::pow(average_data.area.avg(), 2)}}; + } + json_ouput["radius"] = probe_radius; +} + +void Voronota::_to_disk() { + if (output_stream) { + output_stream->flush(); + } +} + +void Voronota::_sample() { + using voronotalt::SimplePoint; + using namespace ranges::cpp20::views; + + // Convert single `Particle` to Voronota's `SimpleSphere` + auto to_sphere = [&](const Particle& p) -> voronotalt::SimpleSphere { + return {{p.pos.x(), p.pos.y(), p.pos.z()}, 0.5 * p.traits().sigma + probe_radius}; + }; + const auto spheres = spc.activeParticles() | transform(to_sphere) | ranges::to_vector; + + voronotalt::RadicalTessellation::Result result; + + if (use_pbc) { + auto to_point = [](const Point& p) -> SimplePoint { return {p.x(), p.y(), p.z()}; }; + const auto corner = 0.5 * spc.geometry.getLength(); + const std::vector box_corners = {to_point(-corner), to_point(corner)}; + voronotalt::RadicalTessellation::construct_full_tessellation(spheres, box_corners, result); + } else { + voronotalt::RadicalTessellation::construct_full_tessellation(spheres, result); + } + + auto areas = result.cells_summaries | transform([](auto& s) { return s.sas_area; }); + const auto total_area = std::accumulate(areas.begin(), areas.end(), 0.0); + average_data.area.add(total_area); + average_data.area_squared.add(total_area * total_area); + + if (output_stream) { + *output_stream << this->getNumberOfSteps() << " " << total_area << "\n"; + } +} + +} // namespace Faunus::analysis diff --git a/src/voronotalt/LICENSE.txt b/src/voronotalt/LICENSE.txt new file mode 100644 index 000000000..970f237db --- /dev/null +++ b/src/voronotalt/LICENSE.txt @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2023 Kliment Olechnovic + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/src/voronotalt/VERSION.txt b/src/voronotalt/VERSION.txt new file mode 100644 index 000000000..c561876cf --- /dev/null +++ b/src/voronotalt/VERSION.txt @@ -0,0 +1 @@ +Voronota-LT version 0.9.2 diff --git a/src/voronotalt/basic_types_and_functions.h b/src/voronotalt/basic_types_and_functions.h new file mode 100644 index 000000000..e82d33ad6 --- /dev/null +++ b/src/voronotalt/basic_types_and_functions.h @@ -0,0 +1,370 @@ +#ifndef VORONOTALT_BASIC_TYPES_AND_FUNCTIONS_H_ +#define VORONOTALT_BASIC_TYPES_AND_FUNCTIONS_H_ + +#include + +//#define VORONOTALT_FP32 +#define VORONOTALT_UI32 + +#ifdef VORONOTALT_FP32 +#define FLOATCONST( v ) v ## f +#define FLOATEPSILON 0.000001f +#define PIVALUE 3.14159265358979323846f +#else +#define FLOATCONST( v ) v +#define FLOATEPSILON 0.0000000001 +#define PIVALUE 3.14159265358979323846 +#endif + +namespace voronotalt +{ + +#ifdef VORONOTALT_FP32 +typedef float Float; +#else +typedef double Float; +#endif + +#ifdef VORONOTALT_UI32 +typedef unsigned int UnsignedInt; +#else +typedef std::size_t UnsignedInt; +#endif + +struct SimplePoint +{ + Float x; + Float y; + Float z; + + SimplePoint() : x(FLOATCONST(0.0)), y(FLOATCONST(0.0)), z(FLOATCONST(0.0)) + { + } + + SimplePoint(const Float x, const Float y, const Float z) : x(x), y(y), z(z) + { + } +}; + +struct SimpleSphere +{ + SimplePoint p; + Float r; + + SimpleSphere() : r(FLOATCONST(0.0)) + { + } + + SimpleSphere(const SimplePoint& p, const Float r) : p(p), r(r) + { + } +}; + +struct SimpleQuaternion +{ + Float a; + Float b; + Float c; + Float d; + + SimpleQuaternion(const Float a, const Float b, const Float c, const Float d) : a(a), b(b), c(c), d(d) + { + } + + SimpleQuaternion(const Float a, const SimplePoint& p) : a(a), b(p.x), c(p.y), d(p.z) + { + } +}; + +inline bool equal(const Float a, const Float b, const Float e) +{ + return (((a-b)<=e) && ((b-a)<=e)); +} + +inline bool equal(const Float a, const Float b) +{ + return equal(a, b, FLOATEPSILON); +} + +inline bool less(const Float a, const Float b) +{ + return ((a+FLOATEPSILON)b); +} + +inline bool less_or_equal(const Float a, const Float b) +{ + return (less(a, b) || equal(a, b)); +} + +inline bool greater_or_equal(const Float a, const Float b) +{ + return (greater(a, b) || equal(a, b)); +} + +inline void set_close_to_equal(Float& a, const Float b) +{ + if(equal(a, b)) + { + a=b; + } +} + +inline Float squared_distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) +{ + const Float dx=(a.x-b.x); + const Float dy=(a.y-b.y); + const Float dz=(a.z-b.z); + return (dx*dx+dy*dy+dz*dz); +} + +inline Float distance_from_point_to_point(const SimplePoint& a, const SimplePoint& b) +{ + return sqrt(squared_distance_from_point_to_point(a, b)); +} + +inline Float squared_point_module(const SimplePoint& a) +{ + return (a.x*a.x+a.y*a.y+a.z*a.z); +} + +inline Float point_module(const SimplePoint& a) +{ + return sqrt(squared_point_module(a)); +} + +inline Float dot_product(const SimplePoint& a, const SimplePoint& b) +{ + return (a.x*b.x+a.y*b.y+a.z*b.z); +} + +inline SimplePoint cross_product(const SimplePoint& a, const SimplePoint& b) +{ + return SimplePoint(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x); +} + +inline SimplePoint point_and_number_product(const SimplePoint& a, const Float k) +{ + return SimplePoint(a.x*k, a.y*k, a.z*k); +} + +inline SimplePoint unit_point(const SimplePoint& a) +{ + return ((equal(squared_point_module(a), FLOATCONST(1.0))) ? a : point_and_number_product(a, FLOATCONST(1.0)/point_module(a))); +} + +inline SimplePoint sum_of_points(const SimplePoint& a, const SimplePoint& b) +{ + return SimplePoint(a.x+b.x, a.y+b.y, a.z+b.z); +} + +inline SimplePoint sub_of_points(const SimplePoint& a, const SimplePoint& b) +{ + return SimplePoint(a.x-b.x, a.y-b.y, a.z-b.z); +} + +inline Float signed_distance_from_point_to_plane(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) +{ + return dot_product(unit_point(plane_normal), sub_of_points(x, plane_point)); +} + +inline int halfspace_of_point(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& x) +{ + const Float sd=signed_distance_from_point_to_plane(plane_point, plane_normal, x); + if(greater(sd, FLOATCONST(0.0))) + { + return 1; + } + else if(less(sd, FLOATCONST(0.0))) + { + return -1; + } + return 0; +} + +inline SimplePoint intersection_of_plane_and_segment(const SimplePoint& plane_point, const SimplePoint& plane_normal, const SimplePoint& a, const SimplePoint& b) +{ + const Float da=signed_distance_from_point_to_plane(plane_point, plane_normal, a); + const Float db=signed_distance_from_point_to_plane(plane_point, plane_normal, b); + if((da-db)==0) + { + return a; + } + else + { + const Float t=da/(da-db); + return sum_of_points(a, point_and_number_product(sub_of_points(b, a), t)); + } +} + +inline Float triangle_area(const SimplePoint& a, const SimplePoint& b, const SimplePoint& c) +{ + return (point_module(cross_product(sub_of_points(b, a), sub_of_points(c, a)))/FLOATCONST(2.0)); +} + +inline Float min_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b) +{ + Float cos_val=dot_product(unit_point(sub_of_points(a, o)), unit_point(sub_of_points(b, o))); + if(cos_valFLOATCONST(1.0)) + { + cos_val=FLOATCONST(1.0); + } + return std::acos(cos_val); +} + +inline Float directed_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b, const SimplePoint& c) +{ + const Float angle=min_angle(o, a, b); + const SimplePoint n=cross_product(unit_point(sub_of_points(a, o)), unit_point(sub_of_points(b, o))); + if(dot_product(sub_of_points(c, o), n)>=0) + { + return angle; + } + else + { + return (PIVALUE*FLOATCONST(2.0)-angle); + } +} + +SimplePoint any_normal_of_vector(const SimplePoint& a) +{ + SimplePoint b=a; + if(!equal(b.x, FLOATCONST(0.0)) && (!equal(b.y, FLOATCONST(0.0)) || !equal(b.z, FLOATCONST(0.0)))) + { + b.x=FLOATCONST(0.0)-b.x; + return unit_point(cross_product(a, b)); + } + else if(!equal(b.y, FLOATCONST(0.0)) && (!equal(b.x, FLOATCONST(0.0)) || !equal(b.z, FLOATCONST(0.0)))) + { + b.y=FLOATCONST(0.0)-b.y; + return unit_point(cross_product(a, b)); + } + else if(!equal(b.x, FLOATCONST(0.0))) + { + return SimplePoint(FLOATCONST(0.0), FLOATCONST(1.0), FLOATCONST(0.0)); + } + else + { + return SimplePoint(FLOATCONST(1.0), FLOATCONST(0.0), FLOATCONST(0.0)); + } +} + +inline bool sphere_intersects_sphere(const SimpleSphere& a, const SimpleSphere& b) +{ + return less(squared_distance_from_point_to_point(a.p, b.p), (a.r+b.r)*(a.r+b.r)); +} + +inline bool sphere_equals_sphere(const SimpleSphere& a, const SimpleSphere& b) +{ + return (equal(a.r, b.r) && equal(a.p.x, b.p.x) && equal(a.p.y, b.p.y) && equal(a.p.z, b.p.z)); +} + +inline bool sphere_contains_sphere(const SimpleSphere& a, const SimpleSphere& b) +{ + return (greater_or_equal(a.r, b.r) && less_or_equal(squared_distance_from_point_to_point(a.p, b.p), (a.r-b.r)*(a.r-b.r))); +} + +inline Float distance_to_center_of_intersection_circle_of_two_spheres(const SimpleSphere& a, const SimpleSphere& b) +{ + const Float cm=distance_from_point_to_point(a.p, b.p); + const Float cos_g=(a.r*a.r+cm*cm-b.r*b.r)/(2*a.r*cm); + return (a.r*cos_g); +} + +inline SimplePoint center_of_intersection_circle_of_two_spheres(const SimpleSphere& a, const SimpleSphere& b) +{ + const SimplePoint cv=sub_of_points(b.p, a.p); + const Float cm=point_module(cv); + const Float cos_g=(a.r*a.r+cm*cm-b.r*b.r)/(2*a.r*cm); + return sum_of_points(a.p, point_and_number_product(cv, a.r*cos_g/cm)); +} + +inline SimpleSphere intersection_circle_of_two_spheres(const SimpleSphere& a, const SimpleSphere& b) +{ + const SimplePoint cv=sub_of_points(b.p, a.p); + const Float cm=point_module(cv); + const Float cos_g=(a.r*a.r+cm*cm-b.r*b.r)/(2*a.r*cm); + const Float sin_g=std::sqrt(1-cos_g*cos_g); + return SimpleSphere(sum_of_points(a.p, point_and_number_product(cv, a.r*cos_g/cm)), a.r*sin_g); +} + +inline bool project_point_inside_line(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b, SimplePoint& result) +{ + const SimplePoint v=unit_point(sub_of_points(b, a)); + const Float l=dot_product(v, sub_of_points(o, a)); + if(l>FLOATCONST(0.0) && (l*l)<=squared_distance_from_point_to_point(a, b)) + { + result=sum_of_points(a, point_and_number_product(v, l)); + return true; + } + return false; +} + +inline bool intersect_segment_with_circle(const SimpleSphere& circle, const SimplePoint& p_in, const SimplePoint& p_out, SimplePoint& result) +{ + const Float dist_in_to_out=distance_from_point_to_point(p_in, p_out); + if(dist_in_to_out>FLOATCONST(0.0)) + { + const SimplePoint v=point_and_number_product(sub_of_points(p_in, p_out), FLOATCONST(1.0)/dist_in_to_out); + const SimplePoint u=sub_of_points(circle.p, p_out); + const SimplePoint s=sum_of_points(p_out, point_and_number_product(v, dot_product(v, u))); + const Float ll=(circle.r*circle.r)-squared_distance_from_point_to_point(circle.p, s); + if(ll>=FLOATCONST(0.0)) + { + result=sum_of_points(s, point_and_number_product(v, FLOATCONST(0.0)-std::sqrt(ll))); + return true; + } + } + return false; +} + +inline Float min_dihedral_angle(const SimplePoint& o, const SimplePoint& a, const SimplePoint& b1, const SimplePoint& b2) +{ + const SimplePoint oa=unit_point(sub_of_points(a, o)); + const SimplePoint d1=sub_of_points(b1, sum_of_points(o, point_and_number_product(oa, dot_product(oa, sub_of_points(b1, o))))); + const SimplePoint d2=sub_of_points(b2, sum_of_points(o, point_and_number_product(oa, dot_product(oa, sub_of_points(b2, o))))); + const Float cos_val=dot_product(unit_point(d1), unit_point(d2)); + return std::acos(std::max(FLOATCONST(-1.0), std::min(cos_val, FLOATCONST(1.0)))); +} + +inline SimpleQuaternion quaternion_product(const SimpleQuaternion& q1, const SimpleQuaternion& q2) +{ + return SimpleQuaternion( + q1.a*q2.a - q1.b*q2.b - q1.c*q2.c - q1.d*q2.d, + q1.a*q2.b + q1.b*q2.a + q1.c*q2.d - q1.d*q2.c, + q1.a*q2.c - q1.b*q2.d + q1.c*q2.a + q1.d*q2.b, + q1.a*q2.d + q1.b*q2.c - q1.c*q2.b + q1.d*q2.a); +} + +inline SimpleQuaternion inverted_quaternion(const SimpleQuaternion& q) +{ + return SimpleQuaternion(q.a, FLOATCONST(0.0)-q.b, FLOATCONST(0.0)-q.c, FLOATCONST(0.0)-q.d); +} + +inline SimplePoint rotate_point_around_axis(const SimplePoint& axis, const Float angle, const SimplePoint& p) +{ + if(squared_point_module(axis)>0) + { + const Float radians_angle_half=(angle*FLOATCONST(0.5)); + const SimpleQuaternion q1=SimpleQuaternion(std::cos(radians_angle_half), point_and_number_product(unit_point(axis), std::sin(radians_angle_half))); + const SimpleQuaternion q2=SimpleQuaternion(FLOATCONST(0.0), p); + const SimpleQuaternion q3=quaternion_product(quaternion_product(q1, q2), inverted_quaternion(q1)); + return SimplePoint(q3.b, q3.c, q3.d); + } + else + { + return p; + } +} + +} + +#endif /* VORONOTALT_BASIC_TYPES_AND_FUNCTIONS_H_ */ diff --git a/src/voronotalt/conversion_to_input.h b/src/voronotalt/conversion_to_input.h new file mode 100644 index 000000000..473530f6c --- /dev/null +++ b/src/voronotalt/conversion_to_input.h @@ -0,0 +1,70 @@ +#ifndef VORONOTALT_CONVERSION_TO_INPUT_H_ +#define VORONOTALT_CONVERSION_TO_INPUT_H_ + +#include + +#include "basic_types_and_functions.h" + +namespace voronotalt +{ + +template +inline SimpleSphere get_sphere_from_ball(const Ball& ball, const Float probe) +{ + return SimpleSphere(SimplePoint(static_cast(ball.x), static_cast(ball.y), static_cast(ball.z)), static_cast(ball.r)+probe); +} + +template +inline void fill_sphere_from_ball(const Ball& ball, const Float probe, SimpleSphere& sphere) +{ + sphere.p.x=static_cast(ball.x); + sphere.p.y=static_cast(ball.y); + sphere.p.z=static_cast(ball.z); + sphere.r=static_cast(ball.r)+probe; +} + +template +inline std::vector get_spheres_from_balls(const BallsContainer& balls, const Float probe) +{ + std::vector result; + result.reserve(balls.size()); + for(typename BallsContainer::const_iterator it=balls.begin();it!=balls.end();++it) + { + result.push_back(get_sphere_from_ball(*it, probe)); + } + return result; +} + +template +inline void fill_spheres_from_balls(const BallsContainer& balls, const Float probe, std::vector& spheres) +{ + spheres.resize(balls.size()); + std::size_t i=0; + for(typename BallsContainer::const_iterator it=balls.begin();it!=balls.end();++it) + { + fill_sphere_from_ball(*it, probe, spheres[i]); + i++; + } +} + +template +inline SimplePoint get_simple_point_from_point(const Point& point) +{ + return SimplePoint(static_cast(point.x), static_cast(point.y), static_cast(point.z)); +} + +template +inline std::vector get_simple_points_from_points(const PointsContainer& points) +{ + std::vector result; + result.reserve(points.size()); + for(typename PointsContainer::const_iterator it=points.begin();it!=points.end();++it) + { + result.push_back(get_simple_point_from_point(*it)); + } + return result; +} + +} + +#endif /* VORONOTALT_CONVERSION_TO_INPUT_H_ */ diff --git a/src/voronotalt/preparation_for_tessellation.h b/src/voronotalt/preparation_for_tessellation.h new file mode 100644 index 000000000..31a6a964f --- /dev/null +++ b/src/voronotalt/preparation_for_tessellation.h @@ -0,0 +1,193 @@ +#ifndef VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ +#define VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ + +#include + +#include "spheres_searcher.h" +#include "time_recorder.h" + +namespace voronotalt +{ + +class PreparationForTessellation +{ +public: + struct Result + { + std::vector spheres_with_periodic_instances; + std::vector periodic_links_of_spheres; + std::vector all_exclusion_statuses; + std::vector< std::vector > all_colliding_ids; + std::vector< std::pair > relevant_collision_ids; + UnsignedInt total_input_spheres; + UnsignedInt total_spheres; + UnsignedInt total_collisions; + + Result() : total_input_spheres(0), total_spheres(0), total_collisions(0) + { + } + }; + + static void prepare_for_tessellation( + const std::vector& spheres, + const std::vector& grouping_of_spheres, + Result& result, + TimeRecorder& time_recorder) + { + prepare_for_tessellation(spheres, grouping_of_spheres, std::vector(), result, time_recorder); + } + + static void prepare_for_tessellation( + const std::vector& input_spheres, + const std::vector& grouping_of_spheres, + const std::vector& periodic_box_corners, + Result& result, + TimeRecorder& time_recorder) + { + time_recorder.reset(); + + result=Result(); + + result.total_input_spheres=input_spheres.size(); + + if(periodic_box_corners.size()>=2) + { + SimplePoint corner_a=periodic_box_corners[0]; + SimplePoint corner_b=periodic_box_corners[0]; + for(UnsignedInt i=1;i(sx)); + const bool useful_x=((sx<0 && (corner_a.x-m.p.x)<=buffer_distance) || (sx>0 && (m.p.x-corner_b.x)<=buffer_distance)); + for(int sy=-1;sy<=1;sy++) + { + m.p.y=o.p.y+(shift.y*static_cast(sy)); + const bool useful_y=((sy<0 && (corner_a.y-m.p.y)<=buffer_distance) || (sy>0 && (m.p.y-corner_b.y)<=buffer_distance)); + for(int sz=-1;sz<=1;sz++) + { + if(sx!=0 || sy!=0 || sz!=0) + { + m.p.z=o.p.z+(shift.z*static_cast(sz)); + const bool useful_z=((sz<0 && (corner_a.z-m.p.z)<=buffer_distance) || (sz>0 && (m.p.z-corner_b.z)<=buffer_distance)); + if(useful_x || useful_y || useful_z) + { + result.spheres_with_periodic_instances.push_back(m); + result.periodic_links_of_spheres.push_back(i); + } + } + } + } + } + } + } + + const std::vector& spheres=(result.spheres_with_periodic_instances.empty() ? input_spheres : result.spheres_with_periodic_instances); + + const UnsignedInt N=spheres.size(); + result.total_spheres=N; + + SpheresSearcher spheres_searcher(spheres); + + time_recorder.record_elapsed_miliseconds_and_reset("init spheres searcher"); + + result.all_exclusion_statuses.resize(N, 0); + + result.all_colliding_ids.resize(N); + for(UnsignedInt i=0;i colliding_ids; + colliding_ids.reserve(100); + + std::vector< std::pair > distances_of_colliding_ids; + distances_of_colliding_ids.reserve(100); + + #pragma omp for + for(UnsignedInt i=0;i=grouping_of_spheres.size() || id_b>=grouping_of_spheres.size() || grouping_of_spheres[id_a]!=grouping_of_spheres[id_b]) + { + if(result.periodic_links_of_spheres.empty() || id_a>=result.periodic_links_of_spheres.size() || id_b>=result.periodic_links_of_spheres.size() || result.periodic_links_of_spheres[id_a]==id_a || result.periodic_links_of_spheres[id_b]==id_b) + { + result.relevant_collision_ids.push_back(std::pair(id_a, id_b)); + } + } + } + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("collect relevant collision indices"); + } +}; + +} + +#endif /* VORONOTALT_PREPARATION_FOR_TESSELLATION_H_ */ diff --git a/src/voronotalt/radical_tessellation.h b/src/voronotalt/radical_tessellation.h new file mode 100644 index 000000000..909147d7a --- /dev/null +++ b/src/voronotalt/radical_tessellation.h @@ -0,0 +1,607 @@ +#ifndef VORONOTALT_RADICAL_TESSELLATION_H_ +#define VORONOTALT_RADICAL_TESSELLATION_H_ + +#include + +#include "preparation_for_tessellation.h" +#include "radical_tessellation_contact_construction.h" +#include "time_recorder.h" + +namespace voronotalt +{ + +class RadicalTessellation +{ +public: + struct ContactDescriptorSummary + { + Float area; + Float arc_length; + Float solid_angle_a; + Float solid_angle_b; + Float pyramid_volume_a; + Float pyramid_volume_b; + Float distance; + UnsignedInt flags; + UnsignedInt id_a; + UnsignedInt id_b; + + ContactDescriptorSummary() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + solid_angle_a(FLOATCONST(0.0)), + solid_angle_b(FLOATCONST(0.0)), + pyramid_volume_a(FLOATCONST(0.0)), + pyramid_volume_b(FLOATCONST(0.0)), + distance(FLOATCONST(0.0)), + flags(0), + id_a(0), + id_b(0) + { + } + + void set(const RadicalTessellationContactConstruction::ContactDescriptor& cd) + { + if(cd.area>FLOATCONST(0.0)) + { + id_a=cd.id_a; + id_b=cd.id_b; + area=cd.area; + arc_length=(cd.sum_of_arc_angles*cd.intersection_circle_sphere.r); + solid_angle_a=cd.solid_angle_a; + solid_angle_b=cd.solid_angle_b; + pyramid_volume_a=cd.pyramid_volume_a; + pyramid_volume_b=cd.pyramid_volume_b; + distance=cd.distance; + flags=cd.flags; + } + } + + void ensure_ids_ordered() + { + if(id_a>id_b) + { + std::swap(id_a, id_b); + std::swap(solid_angle_a, solid_angle_b); + std::swap(pyramid_volume_a, pyramid_volume_b); + } + } + }; + + struct CellContactDescriptorsSummary + { + Float area; + Float arc_length; + Float explained_solid_angle_positive; + Float explained_solid_angle_negative; + Float explained_pyramid_volume_positive; + Float explained_pyramid_volume_negative; + Float sas_area; + Float sas_inside_volume; + UnsignedInt id; + UnsignedInt count; + int stage; + + CellContactDescriptorsSummary() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + explained_solid_angle_positive(FLOATCONST(0.0)), + explained_solid_angle_negative(FLOATCONST(0.0)), + explained_pyramid_volume_positive(FLOATCONST(0.0)), + explained_pyramid_volume_negative(FLOATCONST(0.0)), + sas_area(FLOATCONST(0.0)), + sas_inside_volume(FLOATCONST(0.0)), + id(0), + count(0), + stage(0) + { + } + + void add(const ContactDescriptorSummary& cds) + { + if(cds.area>FLOATCONST(0.0) && (cds.id_a==id || cds.id_b==id)) + { + count++; + area+=cds.area; + arc_length+=cds.arc_length; + explained_solid_angle_positive+=std::max(FLOATCONST(0.0), (cds.id_a==id ? cds.solid_angle_a : cds.solid_angle_b)); + explained_solid_angle_negative+=FLOATCONST(0.0)-std::min(FLOATCONST(0.0), (cds.id_a==id ? cds.solid_angle_a : cds.solid_angle_b)); + explained_pyramid_volume_positive+=std::max(FLOATCONST(0.0), (cds.id_a==id ? cds.pyramid_volume_a : cds.pyramid_volume_b)); + explained_pyramid_volume_negative+=FLOATCONST(0.0)-std::min(FLOATCONST(0.0), (cds.id_a==id ? cds.pyramid_volume_a : cds.pyramid_volume_b)); + stage=1; + } + } + + void add(const UnsignedInt new_id, const ContactDescriptorSummary& cds) + { + if(cds.area>FLOATCONST(0.0)) + { + if(stage==0) + { + id=new_id; + } + add(cds); + } + } + + void compute_sas(const Float r) + { + if(stage==1) + { + sas_area=FLOATCONST(0.0); + sas_inside_volume=FLOATCONST(0.0); + if(arc_length>FLOATCONST(0.0) && !equal(explained_solid_angle_positive, explained_solid_angle_negative)) + { + if(explained_solid_angle_positive>explained_solid_angle_negative) + { + sas_area=((FLOATCONST(4.0)*PIVALUE)-std::max(FLOATCONST(0.0), explained_solid_angle_positive-explained_solid_angle_negative))*(r*r); + } + else if(explained_solid_angle_negative>explained_solid_angle_positive) + { + sas_area=(std::max(FLOATCONST(0.0), explained_solid_angle_negative-explained_solid_angle_positive))*(r*r); + } + sas_inside_volume=(sas_area*r/FLOATCONST(3.0))+explained_pyramid_volume_positive-explained_pyramid_volume_negative; + if(sas_inside_volume>(FLOATCONST(4.0)/FLOATCONST(3.0)*PIVALUE*r*r*r)) + { + sas_area=FLOATCONST(0.0); + sas_inside_volume=explained_pyramid_volume_positive-explained_pyramid_volume_negative; + } + } + else + { + sas_inside_volume=explained_pyramid_volume_positive-explained_pyramid_volume_negative; + } + stage=2; + } + } + + void compute_sas_detached(const UnsignedInt new_id, const Float r) + { + if(stage==0) + { + id=new_id; + sas_area=(FLOATCONST(4.0)*PIVALUE)*(r*r); + sas_inside_volume=(sas_area*r/FLOATCONST(3.0)); + stage=2; + } + } + }; + + struct TotalContactDescriptorsSummary + { + Float area; + Float arc_length; + Float distance; + UnsignedInt count; + + TotalContactDescriptorsSummary() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + distance(FLOATCONST(-1.0)), + count(0) + { + } + + void add(const ContactDescriptorSummary& cds) + { + if(cds.area>FLOATCONST(0.0)) + { + count++; + area+=cds.area; + arc_length+=cds.arc_length; + distance=((distance contacts_summaries; + TotalContactDescriptorsSummary total_contacts_summary; + std::vector cells_summaries; + TotalCellContactDescriptorsSummary total_cells_summary; + std::vector contacts_with_redundancy_summaries_in_periodic_box; + std::vector contacts_with_redundancy_canonical_ids_in_periodic_box; + + Result() : total_spheres(0), total_collisions(0), total_relevant_collisions(0) + { + } + }; + + struct ResultGraphics + { + std::vector contacts_graphics; + }; + + struct GroupedResult + { + std::vector grouped_contacts_representative_ids; + std::vector grouped_contacts_summaries; + std::vector grouped_cells_representative_ids; + std::vector grouped_cells_summaries; + }; + + static void construct_full_tessellation( + const std::vector& input_spheres, + Result& result) + { + TimeRecorder time_recorder; + ResultGraphics result_graphics; + construct_full_tessellation(input_spheres, std::vector(), std::vector(), false, true, result, result_graphics, time_recorder); + } + + static void construct_full_tessellation( + const std::vector& input_spheres, + const std::vector& periodic_box_corners, + Result& result) + { + TimeRecorder time_recorder; + ResultGraphics result_graphics; + construct_full_tessellation(input_spheres, std::vector(), periodic_box_corners, false, true, result, result_graphics, time_recorder); + } + + static void construct_full_tessellation( + const std::vector& input_spheres, + const std::vector& grouping_of_spheres, + const std::vector& periodic_box_corners, + const bool with_graphics, + const bool summarize_cells, + Result& result, + ResultGraphics& result_graphics, + TimeRecorder& time_recorder) + { + PreparationForTessellation::Result preparation_result; + PreparationForTessellation::prepare_for_tessellation(input_spheres, grouping_of_spheres, periodic_box_corners, preparation_result, time_recorder); + + const std::vector& spheres=(preparation_result.spheres_with_periodic_instances.empty() ? input_spheres : preparation_result.spheres_with_periodic_instances); + + time_recorder.reset(); + + result=Result(); + result_graphics=ResultGraphics(); + + result.total_spheres=preparation_result.total_input_spheres; + result.total_collisions=preparation_result.total_collisions; + result.total_relevant_collisions=preparation_result.relevant_collision_ids.size(); + + std::vector possible_contacts_summaries(preparation_result.relevant_collision_ids.size()); + + std::vector possible_contacts_graphics; + if(with_graphics) + { + possible_contacts_graphics.resize(possible_contacts_summaries.size()); + } + + time_recorder.record_elapsed_miliseconds_and_reset("allocate possible contact summaries"); + + #pragma omp parallel + { + RadicalTessellationContactConstruction::ContactDescriptor cd; + cd.contour.reserve(12); + + #pragma omp for + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + number_of_valid_contact_summaries++; + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("count valid contact summaries"); + + std::vector ids_of_valid_pairs; + ids_of_valid_pairs.reserve(number_of_valid_contact_summaries); + + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + ids_of_valid_pairs.push_back(i); + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("collect indices of valid contact summaries"); + + result.contacts_summaries.resize(ids_of_valid_pairs.size()); + + #pragma omp parallel + { + #pragma omp for + for(UnsignedInt i=0;i > map_of_spheres_to_boundary_contacts(result.total_spheres); + + for(UnsignedInt i=0;i=result.total_spheres || cds.id_b>=result.total_spheres) && cds.id_a=result.total_spheres || cds.id_b>=result.total_spheres) + && cds.id_a& candidate_ids_a=map_of_spheres_to_boundary_contacts[sphere_id_a]; + const std::vector& candidate_ids_b=map_of_spheres_to_boundary_contacts[sphere_id_b]; + const std::vector& candidate_ids=(candidate_ids_a.size()<=candidate_ids_b.size() ? candidate_ids_a : candidate_ids_b); + UnsignedInt selected_id=result.contacts_summaries.size(); + for(UnsignedInt j=0;j=result.contacts_summaries.size();j++) + { + const UnsignedInt candidate_id=candidate_ids[j]; + const ContactDescriptorSummary& candidate_cds=result.contacts_summaries[candidate_id]; + if(candidate_cds.id_a0) + { + result.contacts_with_redundancy_summaries_in_periodic_box.swap(result.contacts_summaries); + result.contacts_summaries.reserve(result.contacts_with_redundancy_summaries_in_periodic_box.size()+1-count_of_redundant_contacts_in_periodic_box); + for(UnsignedInt i=0;i=result.contacts_with_redundancy_canonical_ids_in_periodic_box.size() || result.contacts_with_redundancy_canonical_ids_in_periodic_box[i]==i) + { + result.contacts_summaries.push_back(cds); + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("reassign ids in contacts at boundaries"); + } + + for(UnsignedInt i=0;i& grouping_of_spheres, + GroupedResult& grouped_result) + { + TimeRecorder time_recorder; + return group_results(full_result, grouping_of_spheres, grouped_result, time_recorder); + } + + static bool group_results( + const Result& full_result, + const std::vector& grouping_of_spheres, + GroupedResult& grouped_result, + TimeRecorder& time_recorder) + { + time_recorder.reset(); + + grouped_result=GroupedResult(); + + if(!grouping_of_spheres.empty() && grouping_of_spheres.size()==full_result.total_spheres) + { + { + grouped_result.grouped_contacts_representative_ids.reserve(full_result.contacts_summaries.size()/5); + grouped_result.grouped_contacts_summaries.reserve(full_result.contacts_summaries.size()/5); + + std::map< std::pair, UnsignedInt > map_of_groups; + + for(UnsignedInt i=0;i group_id(grouping_of_spheres[cds.id_a], grouping_of_spheres[cds.id_b]); + if(group_id.first!=group_id.second) + { + if(group_id.first>group_id.second) + { + std::swap(group_id.first, group_id.second); + } + UnsignedInt group_index=0; + std::map< std::pair, UnsignedInt >::const_iterator it=map_of_groups.find(group_id); + if(it==map_of_groups.end()) + { + group_index=grouped_result.grouped_contacts_summaries.size(); + grouped_result.grouped_contacts_representative_ids.push_back(i); + grouped_result.grouped_contacts_summaries.push_back(TotalContactDescriptorsSummary()); + map_of_groups[group_id]=group_index; + } + else + { + group_index=it->second; + } + grouped_result.grouped_contacts_summaries[group_index].add(cds); + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("grouped contacts summaries"); + } + + if(!full_result.cells_summaries.empty() && full_result.cells_summaries.size()==grouping_of_spheres.size()) + { + grouped_result.grouped_cells_representative_ids.reserve(full_result.cells_summaries.size()/5); + grouped_result.grouped_cells_summaries.reserve(full_result.cells_summaries.size()/5); + + std::map< int, UnsignedInt > map_of_groups; + + for(UnsignedInt i=0;i::const_iterator it=map_of_groups.find(group_id); + if(it==map_of_groups.end()) + { + group_index=grouped_result.grouped_cells_summaries.size(); + grouped_result.grouped_cells_representative_ids.push_back(i); + grouped_result.grouped_cells_summaries.push_back(TotalCellContactDescriptorsSummary()); + map_of_groups[group_id]=group_index; + } + else + { + group_index=it->second; + } + grouped_result.grouped_cells_summaries[group_index].add(ccds); + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("grouped cells summaries"); + } + } + + return (!grouped_result.grouped_contacts_summaries.empty()); + } +}; + +} + +#endif /* VORONOTALT_RADICAL_TESSELLATION_H_ */ diff --git a/src/voronotalt/radical_tessellation_contact_construction.h b/src/voronotalt/radical_tessellation_contact_construction.h new file mode 100644 index 000000000..911829854 --- /dev/null +++ b/src/voronotalt/radical_tessellation_contact_construction.h @@ -0,0 +1,678 @@ +#ifndef VORONOTALT_RADICAL_TESSELLATION_CONTACT_CONSTRUCTION_H_ +#define VORONOTALT_RADICAL_TESSELLATION_CONTACT_CONSTRUCTION_H_ + +#include + +#include "basic_types_and_functions.h" + +namespace voronotalt +{ + +class RadicalTessellationContactConstruction +{ +public: + struct ContourPoint + { + SimplePoint p; + Float angle; + UnsignedInt left_id; + UnsignedInt right_id; + int indicator; + + + ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) : + p(p), + angle(FLOATCONST(0.0)), + left_id(left_id), + right_id(right_id), + indicator(0) + { + } + }; + + typedef std::vector Contour; + + struct ContactDescriptor + { + Contour contour; + SimpleSphere intersection_circle_sphere; + SimplePoint intersection_circle_axis; + SimplePoint contour_barycenter; + Float sum_of_arc_angles; + Float area; + Float solid_angle_a; + Float solid_angle_b; + Float pyramid_volume_a; + Float pyramid_volume_b; + Float distance; + UnsignedInt flags; + UnsignedInt id_a; + UnsignedInt id_b; + + ContactDescriptor() : + sum_of_arc_angles(FLOATCONST(0.0)), + area(FLOATCONST(0.0)), + solid_angle_a(FLOATCONST(0.0)), + solid_angle_b(FLOATCONST(0.0)), + pyramid_volume_a(FLOATCONST(0.0)), + pyramid_volume_b(FLOATCONST(0.0)), + distance(FLOATCONST(0.0)), + flags(0), + id_a(0), + id_b(0) + { + } + + void clear() + { + id_a=0; + id_b=0; + contour.clear(); + sum_of_arc_angles=FLOATCONST(0.0); + area=FLOATCONST(0.0); + solid_angle_a=FLOATCONST(0.0); + solid_angle_b=FLOATCONST(0.0); + pyramid_volume_a=FLOATCONST(0.0); + pyramid_volume_b=FLOATCONST(0.0); + distance=FLOATCONST(0.0); + flags=0; + } + }; + + struct ContactDescriptorGraphics + { + std::vector outer_points; + SimplePoint barycenter; + SimplePoint plane_normal; + + ContactDescriptorGraphics() + { + } + + void clear() + { + outer_points.clear(); + } + }; + + static bool construct_contact_descriptor( + const std::vector& spheres, + const std::vector& spheres_exclusion_statuses, + const UnsignedInt a_id, + const UnsignedInt b_id, + const std::vector& a_neighbor_collisions, + ContactDescriptor& result_contact_descriptor) + { + result_contact_descriptor.clear(); + if(a_idFLOATCONST(0.0)) + { + bool discarded=false; + bool contour_initialized=false; + bool nostop=true; + { + for(UnsignedInt i=0;i=spheres_exclusion_statuses.size() || spheres_exclusion_statuses[neighbor_id]==0)) + { + const SimpleSphere& c=spheres[neighbor_id]; + if(sphere_intersects_sphere(result_contact_descriptor.intersection_circle_sphere, c) && sphere_intersects_sphere(b, c)) + { + if(sphere_contains_sphere(c, a) || sphere_contains_sphere(c, b)) + { + discarded=true; + } + else + { + const SimplePoint neighbor_ac_plane_center=center_of_intersection_circle_of_two_spheres(a, c); + const SimplePoint neighbor_ac_plane_normal=unit_point(sub_of_points(c.p, a.p)); + 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))); + if(std::abs(cos_val)=result_contact_descriptor.intersection_circle_sphere.r) + { + if(halfspace_of_point(neighbor_ac_plane_center, neighbor_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p)>=0) + { + discarded=true; + } + } + else + { + if(!contour_initialized) + { + result_contact_descriptor.intersection_circle_axis=unit_point(sub_of_points(b.p, a.p)); + init_contour_from_base_and_axis(a_id, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.intersection_circle_axis, result_contact_descriptor.contour); + contour_initialized=true; + } + else + { + nostop=(test_if_contour_is_still_cuttable(a.p, neighbor_ac_plane_center, result_contact_descriptor.contour)); + } + if(nostop) + { + mark_and_cut_contour(neighbor_ac_plane_center, neighbor_ac_plane_normal, neighbor_id, result_contact_descriptor.contour); + if(contour_initialized && result_contact_descriptor.contour.empty()) + { + discarded=true; + } + } + } + } + else + { + if(halfspace_of_point(neighbor_ac_plane_center, neighbor_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p)>0) + { + discarded=true; + } + } + } + } + } + } + } + if(!discarded) + { + if(!contour_initialized) + { + result_contact_descriptor.intersection_circle_axis=unit_point(sub_of_points(b.p, a.p)); + result_contact_descriptor.contour_barycenter=result_contact_descriptor.intersection_circle_sphere.p; + result_contact_descriptor.sum_of_arc_angles=(PIVALUE*FLOATCONST(2.0)); + result_contact_descriptor.area=result_contact_descriptor.intersection_circle_sphere.r*result_contact_descriptor.intersection_circle_sphere.r*PIVALUE; + } + else + { + if(!result_contact_descriptor.contour.empty()) + { + 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); + if(!result_contact_descriptor.contour.empty()) + { + result_contact_descriptor.area=calculate_contour_area(result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour, result_contact_descriptor.contour_barycenter); + } + } + } + + if(result_contact_descriptor.area>FLOATCONST(0.0)) + { + result_contact_descriptor.solid_angle_a=calculate_contour_solid_angle(a, b, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour); + result_contact_descriptor.solid_angle_b=calculate_contour_solid_angle(b, a, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.contour); + 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_aFLOATCONST(0.0)); + } + + static bool construct_contact_descriptor_graphics(const ContactDescriptor& contact_descriptor, const Float length_step, ContactDescriptorGraphics& result_contact_descriptor_graphics) + { + result_contact_descriptor_graphics.clear(); + if(contact_descriptor.area>FLOATCONST(0.0)) + { + result_contact_descriptor_graphics.plane_normal=contact_descriptor.intersection_circle_axis; + const Float angle_step=std::max(std::min(length_step/contact_descriptor.intersection_circle_sphere.r, PIVALUE/FLOATCONST(3.0)), PIVALUE/FLOATCONST(36.0)); + if(contact_descriptor.contour.empty()) + { + const SimplePoint first_point=point_and_number_product(any_normal_of_vector(contact_descriptor.intersection_circle_axis), contact_descriptor.intersection_circle_sphere.r); + result_contact_descriptor_graphics.outer_points.reserve(static_cast((PIVALUE*FLOATCONST(2.0))/angle_step)+2); + result_contact_descriptor_graphics.outer_points.push_back(sum_of_points(contact_descriptor.intersection_circle_sphere.p, first_point)); + for(Float rotation_angle=angle_step;rotation_angle<(PIVALUE*FLOATCONST(2.0));rotation_angle+=angle_step) + { + 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))); + } + result_contact_descriptor_graphics.barycenter=contact_descriptor.intersection_circle_sphere.p; + } + else + { + if(contact_descriptor.sum_of_arc_angles>FLOATCONST(0.0)) + { + result_contact_descriptor_graphics.outer_points.reserve(static_cast(contact_descriptor.sum_of_arc_angles/angle_step)+contact_descriptor.contour.size()+4); + } + else + { + result_contact_descriptor_graphics.outer_points.reserve(contact_descriptor.contour.size()); + } + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + if(pr.angle>angle_step) + { + const SimplePoint first_v=sub_of_points(pr.p, contact_descriptor.intersection_circle_sphere.p); + for(Float rotation_angle=angle_step;rotation_angle(result_contact_descriptor_graphics.outer_points.size()); + result_contact_descriptor_graphics.barycenter.y/=static_cast(result_contact_descriptor_graphics.outer_points.size()); + result_contact_descriptor_graphics.barycenter.z/=static_cast(result_contact_descriptor_graphics.outer_points.size()); + } + } + } + return (!result_contact_descriptor_graphics.outer_points.empty()); + } + +private: + static void init_contour_from_base_and_axis( + const UnsignedInt a_id, + const SimpleSphere& base, + const SimplePoint& axis, + Contour& result) + { + const SimplePoint first_point=point_and_number_product(any_normal_of_vector(axis), base.r*FLOATCONST(1.19)); + const Float angle_step=PIVALUE/FLOATCONST(3.0); + result.push_back(ContourPoint(sum_of_points(base.p, first_point), a_id, a_id)); + for(Float rotation_angle=angle_step;rotation_angle<(PIVALUE*FLOATCONST(2.0));rotation_angle+=angle_step) + { + result.push_back(ContourPoint(sum_of_points(base.p, rotate_point_around_axis(axis, rotation_angle, first_point)), a_id, a_id)); + } + } + + static bool mark_and_cut_contour( + const SimplePoint& ac_plane_center, + const SimplePoint& ac_plane_normal, + const UnsignedInt c_id, + Contour& contour) + { + const UnsignedInt outsiders_count=mark_contour(ac_plane_center, ac_plane_normal, c_id, contour); + if(outsiders_count>0) + { + if(outsiders_countp)>=0) + { + it->left_id=c_id; + it->right_id=c_id; + outsiders_count++; + } + } + return outsiders_count; + } + + static void cut_contour( + const SimplePoint& ac_plane_center, + const SimplePoint& ac_plane_normal, + const UnsignedInt c_id, + Contour& contour) + { + if(contour.size()<3) + { + return; + } + + UnsignedInt i_start=0; + while(i_start=contour.size()) + { + return; + } + + UnsignedInt i_end=(contour.size()-1); + while(i_end>0 && !(contour[i_end].left_id==c_id && contour[i_end].right_id==c_id)) + { + i_end--; + } + + if(i_start==0 && i_end==(contour.size()-1)) + { + i_end=0; + while((i_end+1)0 && contour[i_start-1].left_id==c_id && contour[i_start-1].right_id==c_id) + { + i_start--; + } + } + + if(i_start==i_end) + { + contour.insert(contour.begin()+i_start, contour[i_start]); + i_end=i_start+1; + } + else if(i_starti_end) + { + if(i_start+10) + { + contour.erase(contour.begin(), contour.begin()+i_end); + + } + i_start=contour.size()-1; + i_end=0; + } + + { + const UnsignedInt i_left=((i_start)>0 ? (i_start-1) : (contour.size()-1)); + 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); + } + + { + const UnsignedInt i_right=((i_end+1)=dist_threshold; + } + return cuttable; + } + + static bool restrict_contour_to_circle( + const SimpleSphere& ic_sphere, + const SimplePoint& ic_axis, + const UnsignedInt a_id, + Contour& contour, + Float& sum_of_arc_angles) + { + UnsignedInt outsiders_count=0; + for(UnsignedInt i=0;i0) + { + UnsignedInt insertions_count=0; + { + UnsignedInt i=0; + while(i0 ? (i-1) : contour.size()); + } + } + if(contour.size()<2) + { + contour.clear(); + } + { + UnsignedInt i=0; + while(i0); + } + + static Float calculate_contour_area(const SimpleSphere& ic_sphere, const Contour& contour, SimplePoint& contour_barycenter) + { + Float area=FLOATCONST(0.0); + + contour_barycenter.x=FLOATCONST(0.0); + contour_barycenter.y=FLOATCONST(0.0); + contour_barycenter.z=FLOATCONST(0.0); + for(UnsignedInt i=0;i(contour.size()); + contour_barycenter.y/=static_cast(contour.size()); + contour_barycenter.z/=static_cast(contour.size()); + + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + area+=ic_sphere.r*ic_sphere.r*(pr1.angle-std::sin(pr1.angle))*0.5; + } + } + + return area; + } + + static Float calculate_contour_solid_angle(const SimpleSphere& a, const SimpleSphere& b, const SimpleSphere& ic_sphere, const Contour& contour) + { + Float turn_angle=FLOATCONST(0.0); + + if(!contour.empty()) + { + for(UnsignedInt i=0;i0) ? (i-1) : (contour.size()-1)]; + const ContourPoint& pr1=contour[i]; + const ContourPoint& pr2=contour[(i+1)FLOATCONST(0.0)) + { + SimplePoint d=cross_product(sub_of_points(b.p, a.p), sub_of_points(pr1.p, ic_sphere.p)); + if((pr0.anglePIVALUE && dot_product(d, sub_of_points(pr0.p, pr1.p))>FLOATCONST(0.0))) + { + d=point_and_number_product(d, FLOATCONST(-1.0)); + } + turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, sum_of_points(pr1.p, d), pr2.p)); + } + else if(pr1.angle>FLOATCONST(0.0)) + { + SimplePoint d=cross_product(sub_of_points(b.p, a.p), sub_of_points(pr1.p, ic_sphere.p)); + if((pr1.anglePIVALUE && dot_product(d, sub_of_points(pr2.p, pr1.p))>FLOATCONST(0.0))) + { + d=point_and_number_product(d, FLOATCONST(-1.0)); + } + turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, pr0.p, sum_of_points(pr1.p, d))); + turn_angle+=pr1.angle*(distance_from_point_to_point(a.p, ic_sphere.p)/a.r); + } + else + { + turn_angle+=(PIVALUE-min_dihedral_angle(a.p, pr1.p, pr0.p, pr2.p)); + } + } + } + else + { + turn_angle=(PIVALUE*FLOATCONST(2.0))*(distance_from_point_to_point(a.p, ic_sphere.p)/a.r); + } + + Float solid_angle=((PIVALUE*FLOATCONST(2.0))-turn_angle); + + 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) + +#include "preparation_for_tessellation.h" +#include "simplified_aw_tessellation_contact_construction.h" +#include "time_recorder.h" + +namespace voronotalt +{ + +class SimplifiedAWTessellation +{ +public: + struct ContactDescriptorSummary + { + Float area; + Float arc_length; + Float distance; + UnsignedInt id_a; + UnsignedInt id_b; + + ContactDescriptorSummary() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + distance(FLOATCONST(0.0)), + id_a(0), + id_b(0) + { + } + + void set(const SimplifiedAWTessellationContactConstruction::ContactDescriptor& cd) + { + if(cd.area>FLOATCONST(0.0)) + { + id_a=cd.id_a; + id_b=cd.id_b; + area=cd.area; + arc_length=cd.arc_length; + distance=cd.distance; + } + } + + void ensure_ids_ordered() + { + if(id_a>id_b) + { + std::swap(id_a, id_b); + } + } + }; + + struct TotalContactDescriptorsSummary + { + Float area; + Float arc_length; + Float distance; + UnsignedInt count; + + TotalContactDescriptorsSummary() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + distance(FLOATCONST(-1.0)), + count(0) + { + } + + void add(const ContactDescriptorSummary& cds) + { + if(cds.area>FLOATCONST(0.0)) + { + count++; + area+=cds.area; + arc_length+=cds.arc_length; + distance=((distance contacts_summaries; + TotalContactDescriptorsSummary total_contacts_summary; + + Result() : total_spheres(0), total_collisions(0), total_relevant_collisions(0) + { + } + }; + + struct ResultGraphics + { + std::vector contacts_graphics; + }; + + struct GroupedResult + { + std::vector grouped_contacts_representative_ids; + std::vector grouped_contacts_summaries; + }; + + static void construct_full_tessellation( + const std::vector& spheres, + Result& result) + { + TimeRecorder time_recorder; + ResultGraphics result_graphics; + construct_full_tessellation(spheres, std::vector(), false, result, result_graphics, time_recorder); + } + + static void construct_full_tessellation( + const std::vector& spheres, + const std::vector& grouping_of_spheres, + const bool with_graphics, + Result& result, + ResultGraphics& result_graphics, + TimeRecorder& time_recorder) + { + PreparationForTessellation::Result preparation_result; + PreparationForTessellation::prepare_for_tessellation(spheres, grouping_of_spheres, preparation_result, time_recorder); + + time_recorder.reset(); + + result=Result(); + result_graphics=ResultGraphics(); + + result.total_spheres=preparation_result.total_spheres; + result.total_collisions=preparation_result.total_collisions; + result.total_relevant_collisions=preparation_result.relevant_collision_ids.size(); + + std::vector possible_contacts_summaries(preparation_result.relevant_collision_ids.size()); + + std::vector possible_contacts_graphics; + if(with_graphics) + { + possible_contacts_graphics.resize(possible_contacts_summaries.size()); + } + + time_recorder.record_elapsed_miliseconds_and_reset("allocate possible contact summaries"); + + #pragma omp parallel + { + SimplifiedAWTessellationContactConstruction::ContactDescriptor cd; + cd.neighbor_descriptors.reserve(24); + + #pragma omp for + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + number_of_valid_contact_summaries++; + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("count valid contact summaries"); + + std::vector ids_of_valid_pairs; + ids_of_valid_pairs.reserve(number_of_valid_contact_summaries); + + for(UnsignedInt i=0;iFLOATCONST(0.0)) + { + ids_of_valid_pairs.push_back(i); + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("collect indices of valid contact summaries"); + + result.contacts_summaries.resize(ids_of_valid_pairs.size()); + + #pragma omp parallel + { + #pragma omp for + for(UnsignedInt i=0;i& grouping_of_spheres, + GroupedResult& grouped_result) + { + TimeRecorder time_recorder; + return group_results(full_result, grouping_of_spheres, grouped_result, time_recorder); + } + + static bool group_results( + const Result& full_result, + const std::vector& grouping_of_spheres, + GroupedResult& grouped_result, + TimeRecorder& time_recorder) + { + time_recorder.reset(); + + grouped_result=GroupedResult(); + + if(!grouping_of_spheres.empty() && grouping_of_spheres.size()==full_result.total_spheres) + { + grouped_result.grouped_contacts_representative_ids.reserve(full_result.contacts_summaries.size()/5); + grouped_result.grouped_contacts_summaries.reserve(full_result.contacts_summaries.size()/5); + + std::map< std::pair, UnsignedInt > map_of_groups; + + for(UnsignedInt i=0;i group_id(grouping_of_spheres[cds.id_a], grouping_of_spheres[cds.id_b]); + if(group_id.first!=group_id.second) + { + if(group_id.first>group_id.second) + { + std::swap(group_id.first, group_id.second); + } + UnsignedInt group_index=0; + std::map< std::pair, UnsignedInt >::const_iterator it=map_of_groups.find(group_id); + if(it==map_of_groups.end()) + { + group_index=grouped_result.grouped_contacts_summaries.size(); + grouped_result.grouped_contacts_representative_ids.push_back(i); + grouped_result.grouped_contacts_summaries.push_back(TotalContactDescriptorsSummary()); + map_of_groups[group_id]=group_index; + } + else + { + group_index=it->second; + } + grouped_result.grouped_contacts_summaries[group_index].add(cds); + } + } + } + + time_recorder.record_elapsed_miliseconds_and_reset("grouped contacts summaries"); + } + + return (!grouped_result.grouped_contacts_summaries.empty()); + } +}; + +} + +#endif /* VORONOTALT_SIMPLIFIED_AW_TESSELLATION_H_ */ diff --git a/src/voronotalt/simplified_aw_tessellation_contact_construction.h b/src/voronotalt/simplified_aw_tessellation_contact_construction.h new file mode 100644 index 000000000..7c01fd79f --- /dev/null +++ b/src/voronotalt/simplified_aw_tessellation_contact_construction.h @@ -0,0 +1,782 @@ +#ifndef VORONOTALT_SIMPLIFIED_AW_TESSELLATION_CONTACT_CONSTRUCTION_H_ +#define VORONOTALT_SIMPLIFIED_AW_TESSELLATION_CONTACT_CONSTRUCTION_H_ + +#include +#include +#include + +#include "basic_types_and_functions.h" + +namespace voronotalt +{ + +class SimplifiedAWTessellationContactConstruction +{ +public: + struct ContourPoint + { + SimplePoint p; + UnsignedInt left_id; + UnsignedInt right_id; + + ContourPoint(const SimplePoint& p, const UnsignedInt left_id, const UnsignedInt right_id) : p(p), left_id(left_id), right_id(right_id) + { + } + }; + + typedef std::list Contour; + + struct NeighborDescriptor + { + Float sort_value; + UnsignedInt neighbor_id; + + NeighborDescriptor() : sort_value(FLOATCONST(0.0)), neighbor_id(0) + { + } + + bool operator<(const NeighborDescriptor& d) const + { + return (sort_value outer_points; + SimplePoint barycenter; + }; + + struct ContactDescriptorGraphics + { + std::vector contours_graphics; + + ContactDescriptorGraphics() + { + } + + void clear() + { + contours_graphics.clear(); + } + }; + + struct ContactDescriptor + { + std::list contours; + std::vector neighbor_descriptors; + ContactDescriptorGraphics graphics; + SimpleSphere intersection_circle_sphere; + SimplePoint intersection_circle_axis; + Float area; + Float arc_length; + Float distance; + UnsignedInt id_a; + UnsignedInt id_b; + + ContactDescriptor() : + area(FLOATCONST(0.0)), + arc_length(FLOATCONST(0.0)), + distance(FLOATCONST(0.0)), + id_a(0), + id_b(0) + { + } + + void clear() + { + id_a=0; + id_b=0; + neighbor_descriptors.clear(); + contours.clear(); + graphics.clear(); + area=FLOATCONST(0.0); + arc_length=FLOATCONST(0.0); + distance=FLOATCONST(0.0); + } + }; + + static bool construct_contact_descriptor( + const std::vector& spheres, + const std::vector& spheres_exclusion_statuses, + const UnsignedInt a_id, + const UnsignedInt b_id, + const std::vector& a_neighbor_collisions, + const Float step, + const int projections, + const bool record_graphics, + ContactDescriptor& result_contact_descriptor) + { + result_contact_descriptor.clear(); + if(a_idFLOATCONST(0.0)) + { + bool discarded=false; + { + for(UnsignedInt i=0;i=spheres_exclusion_statuses.size() || spheres_exclusion_statuses[neighbor_id]==0)) + { + const SimpleSphere& c=spheres[neighbor_id]; + if(sphere_intersects_sphere(result_contact_descriptor.intersection_circle_sphere, c) && sphere_intersects_sphere(b, c)) + { + if(sphere_contains_sphere(c, a) || sphere_contains_sphere(c, b)) + { + discarded=true; + } + else + { + const SimplePoint nd_ac_plane_center=center_of_intersection_circle_of_two_spheres(a, c); + const Float dist_from_a_center_to_c_center=distance_from_point_to_point(a.p, c.p); + const SimplePoint nd_ac_plane_normal=point_and_number_product(sub_of_points(c.p, a.p), FLOATCONST(1.0)/dist_from_a_center_to_c_center); + const Float dist_from_a_center_to_nd_ac_plane_center=distance_from_point_to_point(a.p, nd_ac_plane_center); + const Float cos_val=dot_product(unit_point(sub_of_points(result_contact_descriptor.intersection_circle_sphere.p, a.p)), point_and_number_product(sub_of_points(nd_ac_plane_center, a.p), FLOATCONST(1.0)/dist_from_a_center_to_nd_ac_plane_center)); + const int hsi=halfspace_of_point(nd_ac_plane_center, nd_ac_plane_normal, result_contact_descriptor.intersection_circle_sphere.p); + if(std::abs(cos_val)=(result_contact_descriptor.intersection_circle_sphere.r+xl_buffer)) + { + if(hsi>=0) + { + discarded=true; + } + } + else + { + NeighborDescriptor nd; + nd.neighbor_id=neighbor_id; + nd.sort_value=(hsi>0 ? (FLOATCONST(0.0)-xl) : xl); + result_contact_descriptor.neighbor_descriptors.push_back(nd); + } + } + else + { + if(hsi>0) + { + discarded=true; + } + } + } + } + } + } + } + if(!discarded) + { + result_contact_descriptor.intersection_circle_axis=unit_point(sub_of_points(b.p, a.p)); + + { + Contour initial_contour; + init_contour_from_base_and_axis(a_id, result_contact_descriptor.intersection_circle_sphere, result_contact_descriptor.intersection_circle_axis, step, initial_contour); + if(!initial_contour.empty()) + { + result_contact_descriptor.contours.push_back(initial_contour); + } + } + + if(!result_contact_descriptor.neighbor_descriptors.empty()) + { + if(!result_contact_descriptor.contours.empty()) + { + std::sort(result_contact_descriptor.neighbor_descriptors.begin(), result_contact_descriptor.neighbor_descriptors.end()); + for(UnsignedInt i=0;i::iterator jt=result_contact_descriptor.contours.begin(); + while(jt!=result_contact_descriptor.contours.end()) + { + Contour& contour=(*jt); + std::list segments; + if(cut_and_split_contour(a, c, c_id, contour, segments)) + { + if(!contour.empty()) + { + mend_contour(a, b, c, c_id, step, projections, contour); + if(check_contour_intersects_sphere(result_contact_descriptor.intersection_circle_sphere, contour)) + { + ++jt; + } + else + { + jt=result_contact_descriptor.contours.erase(jt); + } + } + else + { + if(!segments.empty()) + { + for(std::list::iterator st=segments.begin();st!=segments.end();++st) + { + mend_contour(a, b, c, c_id, step, projections, (*st)); + } + filter_contours_intersecting_sphere(result_contact_descriptor.intersection_circle_sphere, segments); + if(!segments.empty()) + { + result_contact_descriptor.contours.splice(jt, segments); + } + } + jt=result_contact_descriptor.contours.erase(jt); + } + } + else + { + ++jt; + } + } + } + } + + if(!result_contact_descriptor.contours.empty()) + { + const Float tolerated_deviation=FLOATCONST(0.5); + bool strangely_extended=false; + for(std::list::const_iterator it=result_contact_descriptor.contours.begin();it!=result_contact_descriptor.contours.end() && !strangely_extended;++it) + { + const Contour& contour=(*it); + for(Contour::const_iterator jt=contour.begin();jt!=contour.end() && !strangely_extended;++jt) + { + strangely_extended=strangely_extended || (distance_from_point_to_point(jt->p, a.p)-a.r>tolerated_deviation); + strangely_extended=strangely_extended || (distance_from_point_to_point(jt->p, b.p)-b.r>tolerated_deviation); + } + } + if(strangely_extended) + { + SimplePoint safe_center=point_and_number_product(sum_of_points(a.p, b.p), FLOATCONST(0.5)); + std::list forcibly_shrunk_result=result_contact_descriptor.contours; + for(std::list::iterator it=forcibly_shrunk_result.begin();it!=forcibly_shrunk_result.end();++it) + { + Contour& contour=(*it); + for(Contour::iterator jt=contour.begin();jt!=contour.end();++jt) + { + if((distance_from_point_to_point(jt->p, a.p)-a.r>tolerated_deviation) || (distance_from_point_to_point(jt->p, b.p)-b.r>tolerated_deviation)) + { + jt->p=sum_of_points(safe_center, point_and_number_product(unit_point(sub_of_points(jt->p, safe_center)), std::min(a.r, b.r))); + } + } + } + result_contact_descriptor.contours.swap(forcibly_shrunk_result); + } + } + } + + if(!result_contact_descriptor.contours.empty()) + { + for(std::list::const_iterator it=result_contact_descriptor.contours.begin();it!=result_contact_descriptor.contours.end();++it) + { + const Contour& contour=(*it); + SimplePoint barycenter; + result_contact_descriptor.area+=calculate_area_from_contour(contour, a, b, barycenter); + result_contact_descriptor.arc_length+=calculate_arc_length_from_contour(a_id, contour); + if(record_graphics && result_contact_descriptor.area>FLOATCONST(0.0)) + { + result_contact_descriptor.graphics.contours_graphics.push_back(ContourGraphics()); + result_contact_descriptor.graphics.contours_graphics.back().barycenter=barycenter; + collect_points_from_contour(contour, result_contact_descriptor.graphics.contours_graphics.back().outer_points); + } + } + } + + result_contact_descriptor.distance=distance_from_point_to_point(a.p, b.p); + } + } + } + } + return (result_contact_descriptor.area>FLOATCONST(0.0)); + } + +private: + class HyperboloidBetweenTwoSpheres + { + public: + static inline SimplePoint project_point_on_hyperboloid(const SimplePoint& p, const SimpleSphere& s1, const SimpleSphere& s2) + { + if(s1.r>s2.r) + { + return project_point_on_hyperboloid(p, s2, s1); + } + else + { + const SimplePoint dv=point_and_number_product(sub_of_points(s1.p, s2.p), FLOATCONST(0.5)); + const SimplePoint dv_unit=unit_point(dv); + const SimplePoint c=sum_of_points(s2.p, dv); + const SimplePoint cp=sub_of_points(p, c); + const Float lz=dot_product(dv_unit, cp); + const Float lx=std::sqrt(std::max(squared_point_module(cp)-(lz*lz), FLOATCONST(0.0))); + const Float z=project_point_on_simple_hyperboloid(lx, 0, point_module(dv), s1.r, s2.r); + return sum_of_points(sub_of_points(p, point_and_number_product(dv_unit, lz)), point_and_number_product(dv_unit, z)); + } + } + + static inline Float intersect_vector_with_hyperboloid(const SimplePoint& a, const SimplePoint& b, const SimpleSphere& s1, const SimpleSphere& s2) + { + if(s1.r>s2.r) + { + return intersect_vector_with_hyperboloid(a, b, s2, s1); + } + else + { + const SimplePoint dv=point_and_number_product(sub_of_points(s1.p, s2.p), FLOATCONST(0.5)); + const SimplePoint dv_unit=unit_point(dv); + const SimplePoint c=sum_of_points(s2.p, dv); + + const SimplePoint ca=sub_of_points(a, c); + const Float maz=dot_product(dv_unit, ca); + const SimplePoint cax=sub_of_points(sub_of_points(a, point_and_number_product(dv_unit, maz)), c); + const SimplePoint cax_unit=unit_point(cax); + const Float max=dot_product(cax_unit, ca); + const Float may=FLOATCONST(0.0); + + const SimplePoint cb=sub_of_points(b, c); + const Float mbz=dot_product(dv_unit, cb); + const Float mbx=dot_product(cax_unit, cb); + const Float mby=std::sqrt(std::max(squared_point_module(cb)-mbz*mbz-mbx*mbx, FLOATCONST(0.0))); + + return intersect_vector_with_simple_hyperboloid(SimplePoint(max, may, maz), SimplePoint(mbx, mby, mbz), point_module(dv), s1.r, s2.r); + } + } + + private: + static inline Float project_point_on_simple_hyperboloid(const Float x, const Float y, const Float d, const Float r1, const Float r2) + { + if(r1>r2) + { + return project_point_on_simple_hyperboloid(x, y, d, r2, r1); + } + else + { + const Float r=r2-r1; + return 2*r*std::sqrt(std::max((0-r*r+4*d*d)*(4*x*x+4*y*y+4*d*d-r*r), FLOATCONST(0.0)))/(0-4*r*r+16*d*d); + } + } + + static inline Float intersect_vector_with_simple_hyperboloid(const SimplePoint& a, const SimplePoint& b, const Float d, const Float r1, const Float r2) + { + if(r1>r2) + { + return intersect_vector_with_simple_hyperboloid(a, b, d, r2, r1); + } + else + { + const Float r=r2-r1; + SimplePoint ab=sub_of_points(b, a); + SimplePoint v=unit_point(ab); + const Float k=(4*r*r/((0-4*r*r+16*d*d)*(0-4*r*r+16*d*d))) * (0-r*r+4*d*d) * 4; + const Float m=(4*d*d-r*r)*k/4; + + const Float x0=a.x; + const Float y0=a.y; + const Float z0=a.z; + const Float vx=v.x; + const Float vy=v.y; + const Float vz=v.z; + + const Float t1 = (std::sqrt((k*vy*vy+k*vx*vx)*z0*z0+(-2*k*vy*vz*y0-2*k*vx*vz*x0)*z0+(k*vz*vz-k*k*vx*vx)*y0*y0+2*k*k*vx*vy*x0*y0+(k*vz*vz-k*k*vy*vy)*x0*x0+m*vz*vz-k*m*vy*vy-k*m*vx*vx)-vz*z0+k*vy*y0+k*vx*x0)/(vz*vz-k*vy*vy-k*vx*vx); + + const Float t2 = -(std::sqrt((k*vy*vy+k*vx*vx)*z0*z0+(-2*k*vy*vz*y0-2*k*vx*vz*x0)*z0+(k*vz*vz-k*k*vx*vx)*y0*y0+2*k*k*vx*vy*x0*y0+(k*vz*vz-k*k*vy*vy)*x0*x0+m*vz*vz-k*m*vy*vy-k*m*vx*vx)+vz*z0-k*vy*y0-k*vx*x0)/(vz*vz-k*vy*vy-k*vx*vx); + + const SimplePoint tp1=sum_of_points(a, point_and_number_product(v, t1)); + const SimplePoint tp2=sum_of_points(a, point_and_number_product(v, t2)); + if(greater(t1, FLOATCONST(0.0)) && less(t1, point_module(ab)) && equal(tp1.z, std::sqrt(k*tp1.x*tp1.x+k*tp1.y*tp1.y+m), FLOATCONST(0.000001))) + { + return t1; + } + else if(greater(t2, FLOATCONST(0.0)) && less(t2, point_module(ab)) && equal(tp2.z, std::sqrt(k*tp2.x*tp2.x+k*tp2.y*tp2.y+m), FLOATCONST(0.000001))) + { + return t2; + } + else + { + return FLOATCONST(0.0); + } + } + } + }; + + static void init_contour_from_base_and_axis( + const UnsignedInt a_id, + const SimpleSphere& base, + const SimplePoint& axis, + const Float length_step, + Contour& result) + { + const Float angle_step=std::max(std::min(length_step/base.r, PIVALUE/FLOATCONST(3.0)), PIVALUE/FLOATCONST(36.0)); + const SimplePoint first_point=point_and_number_product(any_normal_of_vector(axis), base.r); + result.push_back(ContourPoint(sum_of_points(base.p, first_point), a_id, a_id)); + for(Float rotation_angle=angle_step;rotation_angle<(PIVALUE*FLOATCONST(2.0));rotation_angle+=angle_step) + { + result.push_back(ContourPoint(sum_of_points(base.p, rotate_point_around_axis(axis, rotation_angle, first_point)), a_id, a_id)); + } + } + + static bool cut_and_split_contour( + const SimpleSphere& a, + const SimpleSphere& c, + const UnsignedInt c_id, + Contour& contour, + std::list& segments) + { + const UnsignedInt outsiders_count=mark_contour(a, c, c_id, contour); + if(outsiders_count>0) + { + if(outsiders_count cuts; + const int cuts_count=cut_contour(a, c, c_id, contour, cuts); + if(cuts_count>0 && cuts_count%2==0) + { + if(cuts_count>2) + { + order_cuts(cuts); + split_contour(contour, cuts, segments); + } + } + } + else + { + contour.clear(); + } + return true; + } + return false; + } + + static UnsignedInt mark_contour( + const SimpleSphere& a, + const SimpleSphere& c, + const UnsignedInt c_id, + Contour& contour) + { + UnsignedInt outsiders_count=0; + for(Contour::iterator it=contour.begin();it!=contour.end();++it) + { + if((distance_from_point_to_point(it->p, c.p)-c.r)<(distance_from_point_to_point(it->p, a.p)-a.r)) + { + it->left_id=c_id; + it->right_id=c_id; + outsiders_count++; + } + } + return outsiders_count; + } + + static UnsignedInt cut_contour( + const SimpleSphere& a, + const SimpleSphere& c, + const UnsignedInt c_id, + Contour& contour, + std::list& cuts) + { + UnsignedInt cuts_count=0; + Contour::iterator it=contour.begin(); + while(it!=contour.end()) + { + if(it->left_id==c_id && it->right_id==c_id) + { + const Contour::iterator left_it=get_left_iterator(contour, it); + const Contour::iterator right_it=get_right_iterator(contour, it); + + if(left_it->right_id!=c_id) + { + const SimplePoint& p0=it->p; + const SimplePoint& p1=left_it->p; + const Float l=HyperboloidBetweenTwoSpheres::intersect_vector_with_hyperboloid(p0, p1, a, c); + cuts.push_back(contour.insert(it, ContourPoint(sum_of_points(p0, point_and_number_product(unit_point(sub_of_points(p1, p0)), l)), left_it->right_id, it->left_id))); + cuts_count++; + } + + if(right_it->left_id!=c_id) + { + const SimplePoint& p0=it->p; + const SimplePoint& p1=right_it->p; + const Float l=HyperboloidBetweenTwoSpheres::intersect_vector_with_hyperboloid(p0, p1, a, c); + cuts.push_back(contour.insert(right_it, ContourPoint(sum_of_points(p0, point_and_number_product(unit_point(sub_of_points(p1, p0)), l)), it->right_id, right_it->left_id))); + cuts_count++; + } + + it=contour.erase(it); + } + else + { + ++it; + } + } + return cuts_count; + } + + static void order_cuts(std::list& cuts) + { + Float sums[2]={FLOATCONST(0.0), FLOATCONST(0.0)}; + for(int i=0;i<2;i++) + { + if(i==1) + { + shift_list(cuts, false); + } + std::list::const_iterator it=cuts.begin(); + while(it!=cuts.end()) + { + std::list::const_iterator next=it; + ++next; + if(next!=cuts.end()) + { + sums[i]+=distance_from_point_to_point((*it)->p, (*next)->p); + it=next; + ++it; + } + else + { + it=cuts.end(); + } + } + } + if(sums[0]& ordered_cuts, + std::list& segments) + { + UnsignedInt segments_count=0; + std::list::const_iterator it=ordered_cuts.begin(); + while(it!=ordered_cuts.end()) + { + std::list::const_iterator next=it; + ++next; + if(next!=ordered_cuts.end()) + { + if((*next)!=get_right_iterator(contour, (*it))) + { + segments_count++; + Contour segment; + Contour::iterator jt=(*it); + do + { + segment.push_back(*jt); + jt=get_right_iterator(contour, jt); + } + while(jt!=(*next)); + segments.push_back(segment); + } + it=next; + ++it; + } + else + { + it=ordered_cuts.end(); + } + } + if(segments_count>0) + { + contour.clear(); + } + return segments_count; + } + + static void mend_contour( + const SimpleSphere& a, + const SimpleSphere& b, + const SimpleSphere& c, + const UnsignedInt c_id, + const Float step, + const int projections, + Contour& contour) + { + Contour::iterator it=contour.begin(); + while(it!=contour.end()) + { + if(it->left_id!=c_id && it->right_id==c_id) + { + const Contour::iterator jt=get_right_iterator(contour, it); + if(jt->left_id==c_id) + { + SimplePoint& p0=it->p; + SimplePoint& p1=jt->p; + for(int e=0;estep) + { + const int leap_distance=static_cast(std::floor(distance/step+FLOATCONST(0.5))); + const Float leap_size=distance/static_cast(leap_distance); + const SimplePoint direction=unit_point(sub_of_points(p1, p0)); + for(int leap=1;leap(leap)))); + for(int e=0;ep)<=shell.r) + { + return true; + } + } + return false; + } + + static void filter_contours_intersecting_sphere(const SimpleSphere& shell, std::list& contours) + { + std::list::iterator it=contours.begin(); + while(it!=contours.end()) + { + if(!check_contour_intersects_sphere(shell, (*it))) + { + it=contours.erase(it); + } + else + { + ++it; + } + } + } + + template + static Iterator get_left_iterator(List& container, const Iterator& iterator) + { + Iterator left_it=iterator; + if(left_it==container.begin()) + { + left_it=container.end(); + } + --left_it; + return left_it; + } + + template + static Iterator get_right_iterator(List& container, const Iterator& iterator) + { + Iterator right_it=iterator; + ++right_it; + if(right_it==container.end()) + { + right_it=container.begin(); + } + return right_it; + } + + template + static void shift_list(List& list, const bool reverse) + { + if(!reverse) + { + list.push_front(*list.rbegin()); + list.pop_back(); + } + else + { + list.push_back(*list.begin()); + list.pop_front(); + } + } + + static bool collect_points_from_contour(const Contour& contour, std::vector& contour_points) + { + contour_points.clear(); + if(!contour.empty()) + { + contour_points.reserve(contour.size()); + for(Contour::const_iterator jt=contour.begin();jt!=contour.end();++jt) + { + contour_points.push_back(jt->p); + } + } + return (!contour_points.empty()); + } + + static Float calculate_area_from_contour(const Contour& contour, const SimpleSphere& sphere1, const SimpleSphere& sphere2, SimplePoint& contour_barycenter) + { + Float area=FLOATCONST(0.0); + if(!contour.empty()) + { + contour_barycenter.x=FLOATCONST(0.0); + contour_barycenter.y=FLOATCONST(0.0); + contour_barycenter.z=FLOATCONST(0.0); + for(Contour::const_iterator jt=contour.begin();jt!=contour.end();++jt) + { + contour_barycenter.x+=jt->p.x; + contour_barycenter.y+=jt->p.y; + contour_barycenter.z+=jt->p.z; + } + contour_barycenter.x/=static_cast(contour.size()); + contour_barycenter.y/=static_cast(contour.size()); + contour_barycenter.z/=static_cast(contour.size()); + + contour_barycenter=HyperboloidBetweenTwoSpheres::project_point_on_hyperboloid(contour_barycenter, sphere1, sphere2); + + for(Contour::const_iterator jt=contour.begin();jt!=contour.end();++jt) + { + Contour::const_iterator jt_next=jt; + ++jt_next; + if(jt_next==contour.end()) + { + jt_next=contour.begin(); + } + area+=triangle_area(contour_barycenter, jt->p, jt_next->p); + } + } + return area; + } + + static Float calculate_arc_length_from_contour(const UnsignedInt a_id, const Contour& contour) + { + Float arc_length=FLOATCONST(0.0); + for(Contour::const_iterator jt=contour.begin();jt!=contour.end();++jt) + { + Contour::const_iterator jt_next=jt; + ++jt_next; + if(jt_next==contour.end()) + { + jt_next=contour.begin(); + } + if(jt->right_id==a_id && jt_next->left_id==a_id) + { + arc_length+=distance_from_point_to_point(jt->p, jt_next->p); + } + } + return arc_length; + } +}; + +} + +#endif /* VORONOTALT_SIMPLIFIED_AW_TESSELLATION_CONTACT_CONSTRUCTION_H_ */ diff --git a/src/voronotalt/spheres_searcher.h b/src/voronotalt/spheres_searcher.h new file mode 100644 index 000000000..52a73a593 --- /dev/null +++ b/src/voronotalt/spheres_searcher.h @@ -0,0 +1,182 @@ +#ifndef VORONOTALT_SPHERES_SEARCHER_H_ +#define VORONOTALT_SPHERES_SEARCHER_H_ + +#include + +#include "basic_types_and_functions.h" + +namespace voronotalt +{ + +class SpheresSearcher +{ +public: + explicit SpheresSearcher(const std::vector& spheres) : spheres_(spheres), box_size_(FLOATCONST(1.0)) + { + box_size_=calculate_grid_box_size(spheres_, box_size_); + + for(UnsignedInt i=0;i(boxes_.size()); + boxes_.push_back(std::vector(1, i)); + } + else + { + boxes_[box_id].push_back(i); + } + } + } + + static Float calculate_grid_box_size(const std::vector& spheres, const Float min_box_size) + { + Float box_size=min_box_size; + for(UnsignedInt i=0;i& colliding_ids, const bool discard_hidden, int& exclusion_status) const + { + colliding_ids.clear(); + exclusion_status=0; + if(central_id=0) + { + const int box_id=map_of_boxes_[index]; + if(box_id>=0) + { + const std::vector& ids=boxes_[box_id]; + for(UnsignedInt i=0;iid)) + { + colliding_ids.clear(); + exclusion_status=1; + return false; + } + else if(!discard_hidden + || !sphere_contains_sphere(central_sphere, candidate_sphere)) + { + colliding_ids.push_back(id); + } + } + } + } + } + } + } + } + } + return (!colliding_ids.empty()); + } + +private: + struct GridPoint + { + int x; + int y; + int z; + + GridPoint() : x(0), y(0), z(0) + { + } + + GridPoint(const SimpleSphere& s, const Float grid_step) + { + init(s, grid_step); + } + + GridPoint(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) + { + init(s, grid_step, grid_offset); + } + + void init(const SimpleSphere& s, const Float grid_step) + { + x=static_cast(s.p.x/grid_step); + y=static_cast(s.p.y/grid_step); + z=static_cast(s.p.z/grid_step); + } + + void init(const SimpleSphere& s, const Float grid_step, const GridPoint& grid_offset) + { + x=static_cast(s.p.x/grid_step)-grid_offset.x; + y=static_cast(s.p.y/grid_step)-grid_offset.y; + z=static_cast(s.p.z/grid_step)-grid_offset.z; + } + + int index(const GridPoint& grid_size) const + { + return ((x>=0 && y>=0 && z>=0 && x& spheres_; + GridPoint grid_offset_; + GridPoint grid_size_; + std::vector map_of_boxes_; + std::vector< std::vector > boxes_; + Float box_size_; +}; + +} + +#endif /* VORONOTALT_SPHERES_SEARCHER_H_ */ diff --git a/src/voronotalt/time_recorder.h b/src/voronotalt/time_recorder.h new file mode 100644 index 000000000..ac7be6cc7 --- /dev/null +++ b/src/voronotalt/time_recorder.h @@ -0,0 +1,29 @@ +#ifndef VORONOTALT_TIME_RECORDER_H_ +#define VORONOTALT_TIME_RECORDER_H_ + +namespace voronotalt +{ + +class TimeRecorder +{ +public: + TimeRecorder() + { + } + + virtual ~TimeRecorder() + { + } + + virtual void reset() + { + } + + virtual void record_elapsed_miliseconds_and_reset(const char* /*message*/) + { + } +}; + +} + +#endif /* VORONOTALT_TIME_RECORDER_H_ */ diff --git a/src/voronotalt/voronotalt.h b/src/voronotalt/voronotalt.h new file mode 100644 index 000000000..d8d651458 --- /dev/null +++ b/src/voronotalt/voronotalt.h @@ -0,0 +1,8 @@ +#ifndef VORONOTALT_VORONOTALT_H_ +#define VORONOTALT_VORONOTALT_H_ + +#include "conversion_to_input.h" +#include "radical_tessellation.h" +#include "simplified_aw_tessellation.h" + +#endif /* VORONOTALT_VORONOTALT_H_ */