Skip to content

Commit

Permalink
Rename cell_vector to lattice_coordinates
Browse files Browse the repository at this point in the history
  • Loading branch information
@unageek
unageek committed Aug 12, 2024
1 parent 09bbe97 commit bfb0881
Show file tree
Hide file tree
Showing 8 changed files with 310 additions and 311 deletions.
293 changes: 145 additions & 148 deletions include/polatory/isosurface/rmt/lattice.hpp
View file

Large diffs are not rendered by default.

42 changes: 42 additions & 0 deletions include/polatory/isosurface/rmt/lattice_coordinates.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
#pragma once

#include <Eigen/Core>
#include <boost/container_hash/hash.hpp>

namespace polatory::isosurface::rmt {

using lattice_coordinates = Eigen::RowVector3i;

struct lattice_coordinates_hash {
std::size_t operator()(const lattice_coordinates& m) const noexcept {
std::size_t seed{};
boost::hash_combine(seed, m(0));
boost::hash_combine(seed, m(1));
boost::hash_combine(seed, m(2));
return seed;
}
};

struct lattice_coordinates_less {
bool operator()(const lattice_coordinates& lhs, const lattice_coordinates& rhs) const noexcept {
return std::make_tuple(lhs(0), lhs(1), lhs(2)) < std::make_tuple(rhs(0), rhs(1), rhs(2));
}
};

using lattice_coordinates_pair = std::pair<lattice_coordinates, lattice_coordinates>;

struct lattice_coordinates_pair_hash {
std::size_t operator()(const lattice_coordinates_pair& pair) const noexcept {
std::size_t seed{};
boost::hash_combine(seed, lattice_coordinates_hash()(pair.first));
boost::hash_combine(seed, lattice_coordinates_hash()(pair.second));
return seed;
}
};

inline lattice_coordinates_pair make_lattice_coordinates_pair(const lattice_coordinates& lc1,
const lattice_coordinates& lc2) {
return lattice_coordinates_less()(lc1, lc2) ? std::make_pair(lc1, lc2) : std::make_pair(lc2, lc1);
}

} // namespace polatory::isosurface::rmt
36 changes: 18 additions & 18 deletions include/polatory/isosurface/rmt/neighbor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,29 +2,29 @@

#include <array>
#include <polatory/isosurface/rmt/edge.hpp>
#include <polatory/isosurface/rmt/types.hpp>
#include <polatory/isosurface/rmt/lattice_coordinates.hpp>

namespace polatory::isosurface::rmt {

inline const std::array<cell_vector, 14> kNeighborCellVectors{{
cell_vector(1, 0, 0), // 0
cell_vector(1, 0, 1), // 1
cell_vector(0, 0, 1), // 2
cell_vector(1, 1, 0), // 3
cell_vector(1, 1, 1), // 4
cell_vector(0, 1, 1), // 5
cell_vector(0, 1, 0), // 6
cell_vector(-1, 0, 0), // 7
cell_vector(-1, 0, -1), // 8
cell_vector(0, 0, -1), // 9
cell_vector(-1, -1, 0), // A
cell_vector(-1, -1, -1), // B
cell_vector(0, -1, -1), // C
cell_vector(0, -1, 0), // D
inline const std::array<lattice_coordinates, 14> kNeighborLatticeCoordinatesDeltas{{
lattice_coordinates(1, 0, 0), // 0
lattice_coordinates(1, 0, 1), // 1
lattice_coordinates(0, 0, 1), // 2
lattice_coordinates(1, 1, 0), // 3
lattice_coordinates(1, 1, 1), // 4
lattice_coordinates(0, 1, 1), // 5
lattice_coordinates(0, 1, 0), // 6
lattice_coordinates(-1, 0, 0), // 7
lattice_coordinates(-1, 0, -1), // 8
lattice_coordinates(0, 0, -1), // 9
lattice_coordinates(-1, -1, 0), // A
lattice_coordinates(-1, -1, -1), // B
lattice_coordinates(0, -1, -1), // C
lattice_coordinates(0, -1, 0), // D
}};

inline cell_vector neighbor(const cell_vector& cv, edge_index ei) {
return cv + kNeighborCellVectors.at(ei);
inline lattice_coordinates neighbor(const lattice_coordinates& lc, edge_index ei) {
return lc + kNeighborLatticeCoordinatesDeltas.at(ei);
}

} // namespace polatory::isosurface::rmt
8 changes: 4 additions & 4 deletions include/polatory/isosurface/rmt/node_list.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
#pragma once

#include <polatory/isosurface/rmt/lattice_coordinates.hpp>
#include <polatory/isosurface/rmt/node.hpp>
#include <polatory/isosurface/rmt/types.hpp>
#include <unordered_map>

namespace polatory::isosurface::rmt {

class node_list : public std::unordered_map<cell_vector, node, cell_vector_hash> {
class node_list : public std::unordered_map<lattice_coordinates, node, lattice_coordinates_hash> {
public:
node* node_ptr(const cell_vector& cv) {
auto it = find(cv);
node* node_ptr(const lattice_coordinates& lc) {
auto it = find(lc);
return it != end() ? &it->second : nullptr;
}
};
Expand Down
121 changes: 61 additions & 60 deletions include/polatory/isosurface/rmt/primitive_lattice.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,20 +11,21 @@
#include <polatory/geometry/bbox3d.hpp>
#include <polatory/geometry/point3d.hpp>
#include <polatory/isosurface/rmt/edge.hpp>
#include <polatory/isosurface/rmt/lattice_coordinates.hpp>
#include <polatory/isosurface/rmt/neighbor.hpp>
#include <polatory/isosurface/rmt/types.hpp>

namespace polatory::isosurface::rmt {

inline constexpr double inv_sqrt2 = 0.5 * std::numbers::sqrt2;

// Primitive vectors of the body-centered cubic lattice.
inline const std::array<geometry::vector3d, 3> kLatticeVectors{
// A basis for the body-centered cubic lattice.
inline const std::array<geometry::vector3d, 3> kLatticeBasis{
geometry::vector3d{inv_sqrt2, 1.0, 0.0}, geometry::vector3d{-inv_sqrt2, 0.0, 1.0},
geometry::vector3d{inv_sqrt2, -1.0, 0.0}};

// Reciprocal primitive vectors of the body-centered cubic lattice.
inline const std::array<geometry::vector3d, 3> kDualLatticeVectors{
// The dual of the basis above, i.e., kDualLatticeBasis.at(i).dot(kLatticeBasis.at(j)) is
// 1.0 if i == j, otherwise 0.0.
inline const std::array<geometry::vector3d, 3> kDualLatticeBasis{
geometry::vector3d{inv_sqrt2, 0.5, 0.5}, geometry::vector3d{0.0, 0.0, 1.0},
geometry::vector3d{inv_sqrt2, -0.5, 0.5}};

Expand All @@ -36,25 +37,25 @@ class primitive_lattice {
resolution_(resolution),
inv_aniso_(aniso.inverse()),
trans_aniso_(aniso.transpose()),
a0_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeVectors[0])),
a1_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeVectors[1])),
a2_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeVectors[2])),
b0_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeVectors[0] / resolution_)),
b1_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeVectors[1] / resolution_)),
b2_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeVectors[2] / resolution_)),
cv_offset_(compute_cv_offset()),
a0_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeBasis[0])),
a1_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeBasis[1])),
a2_(geometry::transform_vector<3>(inv_aniso_, resolution_ * kLatticeBasis[2])),
b0_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeBasis[0] / resolution_)),
b1_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeBasis[1] / resolution_)),
b2_(geometry::transform_vector<3>(trans_aniso_, kDualLatticeBasis[2] / resolution_)),
lc_origin_(compute_lattice_coordinates_origin()),
first_ext_bbox_(compute_extended_bbox(1)),
second_ext_bbox_(compute_extended_bbox(2)) {}

const geometry::bbox3d& bbox() const { return bbox_; }

geometry::point3d cell_node_point(const cell_vector& cv) const {
return (cv_offset_(0) + cv(0)) * a0_ + (cv_offset_(1) + cv(1)) * a1_ +
(cv_offset_(2) + cv(2)) * a2_;
geometry::point3d position(const lattice_coordinates& lc) const {
return (lc_origin_(0) + lc(0)) * a0_ + (lc_origin_(1) + lc(1)) * a1_ +
(lc_origin_(2) + lc(2)) * a2_;
}

std::pair<int, int> first_cell_vector_range(int m1, int m2) const {
geometry::point3d point = (cv_offset_(1) + m1) * a1_ + (cv_offset_(2) + m2) * a2_;
std::pair<int, int> first_lattice_coordinate_range(int lc1, int lc2) const {
geometry::point3d point = (lc_origin_(1) + lc1) * a1_ + (lc_origin_(2) + lc2) * a2_;
geometry::point3d direction = a0_;
const auto& bbox = second_extended_bbox();

Expand All @@ -76,16 +77,16 @@ class primitive_lattice {
max_t = std::min(max_t, t1);
}

auto min = cell_vector_from_point_unrounded(point + min_t * direction)(0);
auto max = cell_vector_from_point_unrounded(point + max_t * direction)(0);
auto min = lattice_coordinates_unrounded(point + min_t * direction)(0);
auto max = lattice_coordinates_unrounded(point + max_t * direction)(0);
return {static_cast<int>(std::floor(min)), static_cast<int>(std::ceil(max))};
}

std::pair<int, int> second_cell_vector_range(int m2) const {
std::pair<int, int> second_lattice_coordinate_range(int lc2) const {
std::vector<geometry::vector3d> vertices;

geometry::vector3d normal = a0_.cross(a1_);
auto d = -normal.dot((cv_offset_(2) + m2) * a2_);
auto d = -normal.dot((lc_origin_(2) + lc2) * a2_);
auto bbox_vertices = second_extended_bbox().corners();
for (auto i = 0; i < 7; i++) {
for (auto j = i + 1; j < 8; j++) {
Expand All @@ -99,80 +100,80 @@ class primitive_lattice {
geometry::vector3d pq = q - p;
auto t = -(normal.dot(p) + d) / normal.dot(pq);
if (t >= 0 && t <= 1) {
vertices.push_back(p + t * pq);
vertices.emplace_back(p + t * pq);
}
}
}

auto min = std::numeric_limits<double>::infinity();
auto max = -std::numeric_limits<double>::infinity();
for (auto v : vertices) {
auto cv = cell_vector_from_point_unrounded(v);
min = std::min(min, cv(1));
max = std::max(max, cv(1));
for (const auto& v : vertices) {
auto lc = lattice_coordinates_unrounded(v);
min = std::min(min, lc(1));
max = std::max(max, lc(1));
}
return {static_cast<int>(std::floor(min)), static_cast<int>(std::ceil(max))};
}

std::pair<int, int> third_cell_vector_range() const {
std::pair<int, int> third_lattice_coordinate_range() const {
auto vertices = second_extended_bbox().corners();
auto min = std::numeric_limits<double>::infinity();
auto max = -std::numeric_limits<double>::infinity();
for (auto v : vertices.rowwise()) {
auto cv = cell_vector_from_point_unrounded(v);
min = std::min(min, cv(2));
max = std::max(max, cv(2));
auto lc = lattice_coordinates_unrounded(v);
min = std::min(min, lc(2));
max = std::max(max, lc(2));
}
return {static_cast<int>(std::floor(min)), static_cast<int>(std::ceil(max))};
}

cell_vector closest_cell_vector(const geometry::point3d& p) const {
return cell_vector_from_point_unrounded(p).array().round().cast<int>();
lattice_coordinates lattice_coordinates_rounded(const geometry::point3d& p) const {
return lattice_coordinates_unrounded(p).array().round().cast<int>();
}

// Returns a positively-oriented tetrahedron containing the given point.
cell_vectors tetrahedron(const geometry::point3d& p) const {
cell_vectors cvs(4, 3);
auto cvd = cell_vector_from_point_unrounded(p);
geometry::vector3d cvd_int = cvd.array().floor();
geometry::vector3d cvd_frac = cvd - cvd_int;
cell_vector cv = cvd_int.cast<int>();
Eigen::Matrix<int, 4, 3, Eigen::RowMajor> tetrahedron(const geometry::point3d& p) const {
auto lcd = lattice_coordinates_unrounded(p);
geometry::vector3d lcd_int = lcd.array().floor();
geometry::vector3d lcd_frac = lcd - lcd_int;
lattice_coordinates lc = lcd_int.cast<int>();
std::array<int, 3> rank = {0, 1, 2};
std::sort(rank.begin(), rank.end(), [&](auto i, auto j) { return cvd_frac(i) < cvd_frac(j); });
std::sort(rank.begin(), rank.end(), [&](auto i, auto j) { return lcd_frac(i) < lcd_frac(j); });

cvs.row(0) = cv;
Eigen::Matrix<int, 4, 3, Eigen::RowMajor> tet;
tet.row(0) = lc;
auto make_class = [](int i, int j, int k) constexpr -> int { return 9 * i + 3 * j + k; };
switch (make_class(rank.at(0), rank.at(1), rank.at(2))) {
case make_class(2, 1, 0):
cvs.row(1) = neighbor(cv, edge::k0);
cvs.row(2) = neighbor(cv, edge::k3);
tet.row(1) = neighbor(lc, edge::k0);
tet.row(2) = neighbor(lc, edge::k3);
break;
case make_class(2, 0, 1):
cvs.row(1) = neighbor(cv, edge::k3);
cvs.row(2) = neighbor(cv, edge::k6);
tet.row(1) = neighbor(lc, edge::k3);
tet.row(2) = neighbor(lc, edge::k6);
break;
case make_class(0, 2, 1):
cvs.row(1) = neighbor(cv, edge::k6);
cvs.row(2) = neighbor(cv, edge::k5);
tet.row(1) = neighbor(lc, edge::k6);
tet.row(2) = neighbor(lc, edge::k5);
break;
case make_class(0, 1, 2):
cvs.row(1) = neighbor(cv, edge::k5);
cvs.row(2) = neighbor(cv, edge::k2);
tet.row(1) = neighbor(lc, edge::k5);
tet.row(2) = neighbor(lc, edge::k2);
break;
case make_class(1, 0, 2):
cvs.row(1) = neighbor(cv, edge::k2);
cvs.row(2) = neighbor(cv, edge::k1);
tet.row(1) = neighbor(lc, edge::k2);
tet.row(2) = neighbor(lc, edge::k1);
break;
case make_class(1, 2, 0):
cvs.row(1) = neighbor(cv, edge::k1);
cvs.row(2) = neighbor(cv, edge::k0);
tet.row(1) = neighbor(lc, edge::k1);
tet.row(2) = neighbor(lc, edge::k0);
break;
default:
POLATORY_UNREACHABLE();
break;
}
cvs.row(3) = neighbor(cv, edge::k4);
return cvs;
tet.row(3) = neighbor(lc, edge::k4);
return tet;
}

// Returns the bounding box that contains all nodes to be clustered.
Expand All @@ -184,20 +185,20 @@ class primitive_lattice {
const geometry::bbox3d& second_extended_bbox() const { return second_ext_bbox_; }

private:
geometry::vector3d cell_vector_from_point_unrounded(const geometry::point3d& p) const {
return {p.dot(b0_) - cv_offset_(0), p.dot(b1_) - cv_offset_(1), p.dot(b2_) - cv_offset_(2)};
geometry::vector3d lattice_coordinates_unrounded(const geometry::point3d& p) const {
return {p.dot(b0_) - lc_origin_(0), p.dot(b1_) - lc_origin_(1), p.dot(b2_) - lc_origin_(2)};
}

geometry::vector3d compute_cv_offset() const {
geometry::vector3d compute_lattice_coordinates_origin() const {
geometry::point3d center = bbox_.center();
return {std::round(center.dot(b0_)), std::round(center.dot(b1_)), std::round(center.dot(b2_))};
}

geometry::bbox3d compute_extended_bbox(int extension) const {
geometry::vectors3d neigh(14, 3);
for (edge_index ei = 0; ei < 14; ei++) {
auto cv = kNeighborCellVectors.at(ei);
neigh.row(ei) = cv(0) * a0_ + cv(1) * a1_ + cv(2) * a2_;
auto lc = kNeighborLatticeCoordinatesDeltas.at(ei);
neigh.row(ei) = lc(0) * a0_ + lc(1) * a1_ + lc(2) * a2_;
}

geometry::vector3d min_ext = extension * 1.01 * neigh.colwise().minCoeff();
Expand All @@ -216,7 +217,7 @@ class primitive_lattice {
const geometry::vector3d b0_;
const geometry::vector3d b1_;
const geometry::vector3d b2_;
const geometry::vector3d cv_offset_;
const geometry::vector3d lc_origin_;
const geometry::bbox3d first_ext_bbox_;
const geometry::bbox3d second_ext_bbox_;
};
Expand Down
Loading

0 comments on commit bfb0881

Please sign in to comment.