From c7178b70efb60dbcaac3c989bc95843f7f906bee Mon Sep 17 00:00:00 2001 From: qingfengxia Date: Wed, 21 Oct 2020 21:55:51 +0100 Subject: [PATCH] fix geometry hash when value is less than 1 (ppp issue 20) --- src/UniqueId/UniqueId.h | 269 ++++++++++++++++++++++------------------ src/brep_faceter.cc | 9 +- 2 files changed, 153 insertions(+), 125 deletions(-) diff --git a/src/UniqueId/UniqueId.h b/src/UniqueId/UniqueId.h index fc9bfe4..0bd20d9 100644 --- a/src/UniqueId/UniqueId.h +++ b/src/UniqueId/UniqueId.h @@ -1,124 +1,147 @@ -#pragma once - -#include -#include -#include - -#define USE_HALF_FLOAT 1 -#if USE_HALF_FLOAT -#include "UniqueId/half.hpp" -#else -#error "in-house double to unit16_t mapping is not yet implemented" -#endif - -/* * - * these functions are used "GeometryPropertyBuilder" to calc shape unique ID - * from 4 key values: volume (for solid) or surface area, plus centre of mass cooordiante (x, y, z) - * all geometry (shapes) are saved as a compound of solids, for each solid readback, - * using methods in "GeometryPropertyBuilder" to calc 4 key values. then calc unique ID - * using this unique ID, to match and locate meta information (json file) for each solid - * - * dependencies: half.hpp, json.hpp, OCCT - * */ - -namespace PPP -{ - namespace Utilities - { - // todo: typedef uint64_t UniqueIdType; - - /// small endianness is more popular in CPU worlds - bool isBigEndian(void) - { - union { - uint32_t i; - char c[4]; - } bint = {0x01020304}; - - return bint.c[0] == 1; - } - - /// for unsigned int only - template static T round(T value, T tol) - { - T remained = value % (tol); - T r = value / (tol); - if (remained >= (tol / 2)) - r += 1; - return r * (tol); - } - - /// map from double float to half precision float (binary16) - /// with EPS approximately 0.001, so this mapping uses log10 not equal space mapping - /// bit_mask can be used to further reduced EPS - /// then from half float to underneath uint16_t by reinterpret_cast - /// this function use third-party lib: HalfFloat, - /// if the value is close to zero (zeroThreshold = 1e-3), then set it as zero, - std::uint16_t double2uint16(double value) - { - static double zeroThreshold = 1e-4; - std::uint16_t round_precision = 0x0008; // mask out extra significant bits - if (std::fabs(value) < zeroThreshold) - value = 0.0; -#if USE_HALF_FLOAT - using namespace half_float; - half hf(static_cast(value)); - std::uint16_t* p = reinterpret_cast(&hf); - std::uint16_t ret = round(*p, round_precision); - return ret; -#endif - } - - - /// due to unlikely floating point error, calculated ID can be different after rounding - std::vector nearbyIds(uint64_t id) - { - std::vector nids; - // todo: - return nids; - } - - // approximate equal by float point data comparison - bool uniqueIdEqual(uint64_t id1, std::uint64_t id2) - { - std::uint16_t round_precision = 0x0008; // defined in `double2uint16()` - for (size_t i = 0; i < 4; i++) - { - uint16_t i1 = (id1 << i * 16) & 0xFFFF; - uint16_t i2 = (id2 << i * 16) & 0xFFFF; - // todo: this does not consider overflow - if (i1 > i2 + round_precision || i1 < i2 - round_precision) - { - return false; - } - } - return true; - } - - /// todo: endianess problem. - /// this should generate unique ID for a vector of 4 double values - /// this is not universal unique ID (UUID) yet - std::uint64_t uniqueId(const std::vector values) - { - std::uint64_t ret = 0x00000000; - assert(values.size() == 4); - - for (size_t i = 0; i < values.size(); i++) - { - std::uint64_t tmp = double2uint16(values[i]); - ret ^= (tmp << i * 16); - } - return ret; - } - - /// from float array to int[] then into a hash value - /// using cm3, cm as unit, assuming geometry has mm as unit - std::uint64_t geometryUniqueID(double volume, std::vector centerOfMass) - { - const double scale = 1e2; // should equal to that in GeometryPropertyBuilder class in parallel-preprocessor - double gp = volume / (scale * scale * scale); - std::vector pv{gp, centerOfMass[0] / scale, centerOfMass[1] / scale, centerOfMass[2] / scale}; - return uniqueId(pv); - } - } // namespace Utilities +#pragma once + +#include +#include +#include + +#define USE_HALF_FLOAT 1 +#if USE_HALF_FLOAT +#include "UniqueId/half.hpp" +#else +#error "in-house double to unit16_t mapping is not yet implemented" +#endif + +/* * + * these functions are used "GeometryPropertyBuilder" to calc shape unique ID + * from 4 key values: volume (for solid) or surface area, plus centre of mass cooordiante (x, y, z) + * all geometry (shapes) are saved as a compound of solids, for each solid readback, + * using methods in "GeometryPropertyBuilder" to calc 4 key values. then calc unique ID + * using this unique ID, to match and locate meta information (json file) for each solid + * + * dependencies: half.hpp, json.hpp, OCCT + * */ + +namespace PPP +{ + + namespace Utilities + { + + typedef uint64_t UniqueIdType; + + /// detect byte order, alternatively, GCC has `__BYTE_ORDER__` macro + /// little endianness is more popular in CPU worlds +#ifdef __clang__ + // constexpr for this function is not supported by CLANG version 6, 10 + // `error: constexpr function never produces a constant expression` + bool isBigEndian(void) +#else + constexpr bool isBigEndian(void) +#endif + { + union + { + uint32_t i; + char c[4]; + } bint = {0x01020304}; + + return bint.c[0] == 1; + } + + /// for unsigned int only + template static T round(T value, T tol) + { + T remained = value % (tol); + T r = value / (tol); + if (remained >= (tol / 2)) + r += 1; + return r * (tol); + } + + /** + * because half precision max number is 65000 + * CAD ususally use mm as length unit, so volume would overflow for half precision + * scale to deci-meter, or metre, kilometer, depends on geometry size + * */ + const double LENGTH_SCALE = 0.01; + // if the value is close to zero after scaling, regarded as zero in comparison + const static double ZERO_THRESHOLD = 1; + /// rounding: mask out east significant bits for approximately comparison + const std::uint16_t ROUND_PRECISION_MASK = 0x0008; + + + /// map from double float to half precision float (binary16) + /// with EPS approximately 0.001 (when the value is near 1), the max value 65000 + /// bit_mask can be used to further reduced EPS + /// IEEE754 the least signicant bit at the lower bit when in register + /// then from half float to underneath uint16_t by `reinterpret_cast` + /// this function use third-party lib: HalfFloat, + /// Note: if the value is close to zero (LENGTH_ZERO_THRESHOLD), then set it as zero, + std::uint16_t double2uint16(double value) + { + std::uint16_t round_precision = ROUND_PRECISION_MASK; + if (std::fabs(value) < ZERO_THRESHOLD) + value = 0.0; +#if USE_HALF_FLOAT + using namespace half_float; + half hf(static_cast(value)); + std::uint16_t* p = reinterpret_cast(&hf); + std::uint16_t ret = round(*p, round_precision); + return ret; +#endif + } + + + /// due to unlikely floating point error, calculated Id can be different after rounding + std::vector nearbyIds(UniqueIdType) + { + std::vector nids; + throw "todo: not completed"; + return nids; + } + + /// approximate equal by float point data comparison + /// endianess neutral + bool uniqueIdEqual(uint64_t id1, std::uint64_t id2) + { + for (size_t i = 0; i < 4; i++) + { + uint16_t i1 = (id1 << i * 16) & 0xFFFF; + uint16_t i2 = (id2 << i * 16) & 0xFFFF; + // todo: this does not consider overflow + // NaN half precision should works + if (i1 > i2 + ROUND_PRECISION_MASK || i1 < i2 - ROUND_PRECISION_MASK) + { + return false; + } + } + return true; + } + + /// used by GeometryPropertyBuilder class in parallel-preprocessor + /// assuming native endianess, always calculate Id using the same endianness + /// this should generate unique Id for a vector of 4 double values + /// this is not universal unique Id (UUID) yet, usurally std::byte[16] + UniqueIdType uniqueId(const std::vector values) + { + std::uint64_t ret = 0x00000000; + assert(values.size() == 4); + + for (size_t i = 0; i < values.size(); i++) + { + std::uint64_t tmp = double2uint16(values[i]); + ret ^= (tmp << i * 16); + } + return ret; + } + + /// from float array to int[] then into a hash value + UniqueIdType geometryUniqueId(double volume, std::vector centerOfMass) + { + double gp = volume * (LENGTH_SCALE * LENGTH_SCALE * LENGTH_SCALE); + std::vector pv{gp, centerOfMass[0] * LENGTH_SCALE, centerOfMass[1] * LENGTH_SCALE, + centerOfMass[2] * LENGTH_SCALE}; + return uniqueId(pv); + } + } // namespace Utilities } // namespace PPP \ No newline at end of file diff --git a/src/brep_faceter.cc b/src/brep_faceter.cc index 949f833..0cf3d03 100644 --- a/src/brep_faceter.cc +++ b/src/brep_faceter.cc @@ -140,7 +140,7 @@ void perform_faceting(const TopoDS_Face &face, const FacetingTolerance& facet_to std::uint64_t calculate_unique_id(const TopoDS_Shape &shape) { GProp_GProps v_props; BRepGProp::VolumeProperties(shape, v_props); - return PPP::Utilities::geometryUniqueID(v_props.Mass(), + return PPP::Utilities::geometryUniqueId(v_props.Mass(), {v_props.CentreOfMass().X(), v_props.CentreOfMass().Y(), v_props.CentreOfMass().Z()}); @@ -225,7 +225,12 @@ void facet_all_volumes(const TopTools_HSequenceOfShape &shape_list, if (!material.empty()) { material_volumes[material].push_back(vol); } else { - std::cout << "No material found for ID: " << uniqueID << std::endl; + GProp_GProps v_props; + BRepGProp::VolumeProperties(shape, v_props); + std::cout << "Shape sequence "<< i << " : No material found for ID: " << uniqueID << std::endl; + std::cout << "with volume: " << v_props.Mass() << ", CentreOfMass: " + << v_props.CentreOfMass().X() << ", " << v_props.CentreOfMass().Y() + << ", " << v_props.CentreOfMass().Z() << std::endl; } } }